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FOREWORD 

The  present  final  report  for  AFOSR  contract  number  FA9550-09-C-0135  describes  the  results  of  an 
STTR  Phase  I  collaboration  between  Nielsen  Engineering  and  Research  (NEAR),  Inc.  and  the  Ohio 
State  University  (OSU).  Dr.  Patrick  Reisenthel  (NEAR)  was  the  Principal  Investigator  for  this  effort. 
The  OSU  team  was  led  by  Prof.  Theodore  Allen.  The  present  report  includes  contributions  by  Mr. 
Daniel  Lesieutre  (NEAR)  and  Mr.  Soo  Ho  Lee  (OSU),  and  is  a  comprehensive  stand-alone  report  of 
the  activities  conducted  as  part  of  this  STTR  on  the  development  of  novel  multifidelity  modeling  and 
optimization  approaches  applicable  to  multidisciplinary  design  and  analysis  of  aerospace  vehicles. 

In  this  report,  references  to  “the  Team”  or  “our  Team”  pertain  to  the  joint  NEAR-OSU  partnership.  In 
addition  to  this  report,  and  pending  further  results  to  be  developed  as  part  of  Soo  Ho  Lee's  Ph.D. 
dissertation,  our  Team  intends  to  submit  at  least  two  papers  for  publication.  The  first  of  these  papers, 
targeted  for  the  Journal  of  Global  Optimization,  will  focus  on  the  single-fidelity  version  of  radial  basis 
sequential  optimization.  The  second  paper  will  be  concerned  with  the  multifidelity  version  and  will  be 
submitted  to  Structural  and  Multidisciplinary >  Optimization } 

1.  SUMMARY 

The  development  of  revolutionary  aerospace  vehicles  presents  significant  challenges.  Cross-disciplinary 
integration  offers  potentially  important  advantages  early  in  the  design  cycle,  and  to  do  so  accurately  and 
efficiently  requires  merging  multiple  levels  of  fidelity  within  an  engineering  discipline,  as  well  as 
understanding  the  necessary  degree  of  coupling  between  disciplines. 

To  address  these  demands  our  Team  has  developed  methods  and  code  for  a  radial  basis  function-based 
optimization  alternative  to  multifidelity  Kriging  optimization.  Like  the  Kriging-based  method,  our 
code  involves  more  than  a  single  type  of  simulation  model.  The  approach  can  integrate  data  from 
multiple  disciplines  and  takes  into  account  both  the  location  (in  design  space)  and  the  fidelity  of  further 
data  acquisition/infill. 

Specifically,  a  study  of  radial  basis  functions  in  the  context  of  multifidelity  optimization  is  presented, 
including  comparisons  to  prior  modeling  options  and  alternatives.  We  describe  a  general  multifidelity 
framework  and  specific  methods,  including  hybrid  methods  whose  convergence  is  established. 

Numerical  results  on  test  problems  are  used  to  compare  sequential  radial  basis  optimization  (SRBO) 
versus  other  methods.  The  results  indicate  that  methods  based  on  radial  basis  functions  can  compete  with 
Kriging  methods  in  relation  to  efficiency  and  repeated  identification  of  global  optimum  solutions.  As  a 
result,  our  proposed  SRBO  might  be  preferred  because  of  reduced  computational  overhead,  improved 
numerical  stability,  and  flexibility  in  terms  of  modeling  uncertainty  intervals.  The  example  of  an 
aeroelastic  wing  design  using  two  levels  of  fidelity  and  two  disciplines  is  used  to  demonstrate  the 
feasibility  of  the  proposed  approach.  The  report  concludes  with  a  discussion  of  the  results  and 
opportunities  for  future  research. 

Overall,  the  proposed  methods  offer  substantially  reduced  computational  cost  to  arrive  at  the  optimum 
solution.  As  a  result,  there  is  a  potential  for  engineers  and  researchers  to  solve  more  design  problems  in 
less  time  and  achieve  higher  quality  results  at  reduced  costs. 

2.  INTRODUCTION 

The  development  of  increasingly  more  complex  and  integrated  aerospace  vehicles  that  are  capable  of  a 
variety  of  missions  is  leading  engineers  and  scientists  to  consider  nonconventional  airframes,  new 
structural  and  propulsion  concepts  and  their  interactions  early  in  the  design  cycle.  The  technical 


1  We  also  anticipate  at  least  one  paper  presentation  at  an  AIAA  meeting  on  the  aeroelastic  wing  design  problem. 
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challenges  faced  in  the  design  of  these  revolutionary  air  vehicles  range  from  considerations  of  extreme 
environments  (e.g.,  hypersonic/transatmospheric  vehicles)  to  large  dynamic  shape  variations  of  micro 
UAVs.  This  greater  range  of  design  options  must  be  examined  in  cross-disciplinary  terms,  where 
individual  disciplines  are  allowed  to  challenge  the  constraints  of  other  disciplines.  Cross-disciplinary 
design  space  exploration  is  needed  to  reduce  uncertainty  and  increase  knowledge  in  order  to  make 
better  decisions.  The  large  dimensionality  of  such  a  design  space  traditionally  requires  the  use  of  low- 
fidelity  methods  at  the  conceptual  design  stage,  but  there  is  a  need  to  incorporate  high-fidelity  analysis 
tools  early  in  the  design  process.  Mixed  fidelity  models  are  an  enabling  technology,  allowing  the 
design  of  individual  system  components  to  be  understood  within  the  broader  context  of  a  full  system. 
Increased  computational  power  and  the  proliferation  of  new  unmanned  aerial  vehicle  designs  have 
been  driving  forces  in  the  development  of  variable  fidelity  frameworks  over  the  last  decade.  Y et,  a 
number  of  fundamental  issues  in  multifidelity  computing  remain  unresolved. 

Advanced  aerospace  vehicle  design  methods  must  fuse  data  from  multiple  sources  and  various  levels 
of  fidelity.  In  a  multifidelity  framework,  the  resolution  and  fidelity  of  simulations  should  be  tailored 
to  the  requirements  of  the  analysis.  There  are  many  factors  to  consider,  including  the  spatial 
resolution,  temporal  resolution,  physical  processes  being  modeled,  the  number  of  objects,  number  of 
attributes  of  each  object,  and  the  degree  of  interaction  between  these  objects  (Gardner  and  Hennigan, 
2003). 

High-fidelity  methods  are  typically  more  expensive,  and  must  be  used  sparingly  out  of  practical  necessity. 
Other  factors  to  consider  are  that  high-fidelity  analyses  are  assumed  to  be  more  accurate,  but  may  not 
always  be,  and  that  there  are  often  difficulties  related  to  obtaining  good  gradients.  For  example,  there 
might  be  integer/discrete  variables  associated  with  gross  topological  changes,  e.g.,  number  of  compressor 
stages,  number  of  engines,  control  surfaces  and  so  on.  These  amount  to  major  discontinuities  (“cliff 
edges,”  in  Jarrett' s,  2006,  terminology)  which  computationally  efficient  gradient-based  design 
optimization  methods  cannot  cross. 

An  obvious  advantage  of  low-fidelity  methods  is  that  they  are  computationally  inexpensive  and  fast. 

More  rarely  mentioned  are  the  potential  benefits  of  low- fidelity  methods  in  a  multifidelity  framework, 
namely  that  they  enable  effective  design  exploration,  not  merely  because  of  their  speed,  but  as  an  aid  to 
escape  noisy  local  optima  due  to  the  nonsmooth  nature  of  the  design  space,  a  common  occurrence  in  the 
case  in  high-fidelity  analyses.  Low-fidelity  methods  enable  safe  transit  between  local  optima,  effectively 
tunneling  back  and  forth  through  the  cliff  edges  that  discrete  variable  problems  bring  about  (Jarrett, 

2006).  The  disadvantages  of  low-fidelity  tools  are  that  they  may  produce  inaccurate  results.  Also,  low- 
fidelity  tools  may  not  work  for  unconventional  designs  or  strongly  nonlinear  regimes  (Wintzer  et 
al.,  2006),  or  may  not  be  able  to  produce  the  quantities  of  interest  (e.g.,  aerodynamic  drag). 

A  significant  issue  in  multidisciplinary  frameworks  is  the  difficulty  in  specifying  consistency  between 
models.  A  low-fidelity  and  high-fidelity  model  are  said  to  be  weakly  consistent  if  “the  projection  of  the 
state  of  the  higher  resolution  model  to  the  space  of  the  lower  resolution  model  is  sufficiently  close  to  the 
state  of  the  low  resolution  model”  (Gardner  and  Hennigan,  2003).  In  many  applications,  lower  fidelity 
refers  to  lower  resolution,  but  also  different  methodologies  altogether  involving  significant 
approximations,  not  only  missing  detail  but  also  missing  physics  (e.g.,  absence  of  chemical  reactions, 
aeroelastic/unsteady  effects,  viscosity,  or  turbulence).  In  such  cases,  one  cannot  guarantee  even  weak 
consistency. 

True  integration  of  the  results  from  codes  of  varying  fidelity  poses  a  number  of  challenges.  One  of  them 
is  the  fact  that  prediction  codes  of  varying  fidelity  levels  use  different  input  parameterizations  and, 
therefore,  operate  in  different  spaces  (Robinson  et  al.  2006a, b).  Another  challenge  is  the  development  of 
effective  strategies  when  low-  and  high-fidelity  models  are  not  consistent,  or  when  weak  consistency  is 
satisfied  only  in  a  small  region  of  the  design  space,  therefore  limiting  the  efficiency  of  current 
multifidelity  frameworks  and,  potentially,  the  breadth  of  design  options  being  considered.  The  question 
of  optimal  sampling  of  the  design  space  given  limited  resources  is  a  critical  one  that  is  fundamentally 
linked  to  the  characterization  of  uncertainty  and  uncertainty  requirements.  Specifically,  there  is  a  need  to 
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develop  rigorous  approaches  addressing  not  only  where  to  sample  the  design  space,  but  also  at  what  level 
of  fidelity. 

In  the  past  decade,  there  have  been  numerous  implementations  of  variable  fidelity  ideas  based  on  the 
concept  of  “bridge  functions”  as  a  re-anchoring  framework  to  correct  low-fidelity  analyses  to 
approximate  the  results  of  high-fidelity  analyses.  The  “beta”  factor  initially  proposed  by  Chang  et  al. 
(1993)  is  introduced  as  a  local  multiplicative  factor,  although  many  researchers  (e.g.,  Hatanaka  et 
al.,  2006,  Ghoreyshi  et  al.,  2008)  have  found  it  advantageous  to  implement  the  additive  form  of  the 
correction.  These  schemes  incorporate  Taylor  expansion-based  corrections  between  the  low-  and  high- 
fidelity  models,  requiring  the  design  space  to  be  smooth,  and  require  updates  where  both  low-  and  high- 
fidelity  data  must  be  available  at  the  same  point.  These  factors  tend  to  result  in  relatively  frequent  high- 
fidelity  updates  and,  therefore,  only  modest  improvements  in  computational  savings  over  conventional 
optimization.  While  useful  for  establishing  provable  convergence  to  high-fidelity  results/optimization, 
the  local  correction  method  has  a  tendency  to  limit  the  step  sizes  taken  in  optimization  and,  more 
generally,  curtails  the  size  of  the  design  space  exploration. 

As  pointed  out  in  Forrester  et  al.  (2009),  efficient  global  optimization  revolves  around  being  able  to 
successfully  balance  exploration  and  exploitation.  Thus,  there  is  a  need  for  a  multifidelity  framework 
allowing  bold  optimization  steps.  The  present  work  considers  both  hierarchical2  and  integrated3  variable- 
fidelity  metamodeling  methods  whose  uncertainty  is  reduced  with  each  evaluation.  In  this  work,  a  radial 
basis  function  embodiment  of  these  ideas  is  used  to  benchmark  performance  versus  other  methods  and  to 
demonstrate  efficient  multifidelity  multidisciplinary  optimization. 

3.  OBJECTIVES 

Aerospace  finite  element  method  (FEM)  and  computational  fluid  dynamics  (CFD)  simulations  continue 
to  increase  in  fidelity/realism.  Yet,  even  as  computer  speeds  increase,  the  most  realistic  models  (e.g.,  of 
turbulent  flows)  are  still  relatively  slow  computationally,  making  thorough  optimization  difficult.  In  fact, 
the  gap  between  the  computational  speed  of  low-fidelity  models  and  high-fidelity  models  continues  to 
increase  also.  As  a  result,  methods  for  multifidelity  optimization  which  leverage  efficient  low-fidelity 
models  promise  to  enable  thorough  optimization. 

In  the  years  following  Jones,  Schonlau,  and  Welch  (1998),  there  has  been  growing  interest  in 
optimization  using  global  metamodels,  e.g.,  Kriging  or  radial  basis  functions.  Unlike  other  design  of 
experiments  (DOE)  methods  like  sequential  response  surface  methods,  global  metamodels  can 
integrate  and  generate  predictions  over  the  entire  design  space.  As  a  result,  they  can  integrate  an 
entire  data  set  from  an  entire  design  optimization  process  which  can  involve  human  participation  in 
the  process.  Nonglobal  metamodeling  approaches  for  multifidelity  optimization  are  described  in 
Singh  and  Grandhi  (2010).  Advances  in  global  metamodeling  using  radial  basis  functions  include 
Gutmann  (2001),  Holmstrom  (2008),  Floudas  and  Gounaris  (2009),  and  Antanas  (2010).  Our  Team's 
selection  of  the  radial  basis  function  approach  is  based  on  its  suitability  for  use  in  efficient  global 
optimization  and  its  extension  to  the  multifidelity  context  (see  Section  4). 

The  broad  objective  of  this  work  is  to  pioneer  innovative  methods  for  representing,  managing,  and 
fusing  information  of  various  levels  of  fidelity  within  an  engineering  discipline  and  across  multiple 
disciplines  for  a  wide-range  of  analysis  and  design  tools.  The  specific  objectives  of  this  study  include 
(1)  developing  a  multifidelity  multidisciplinary  framework  through  alternative  formulations  to  previous 
methods,  computational  results,  and  rigorous  mathematical  results,  (2)  extending  Huang  et  al.'s  (2006) 
multifidelity  design  space  sampling  strategy  to  manage  system  updates  within  integrated  and  hierarchical 
global  metamodeling  approaches,  (3)  extending  the  rigorous  convergence  results  from  Schonlau  ( 1 997)  to 
multifidelity  optimization  in  the  context  of  radial  basis  function  methods  and  Kriging  models, 


2  also  referred  to  in  this  work  as  “cascading  model” 

3  referred  to  in  this  work  as  the  “augmented  dimensionality  model” 
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(4)  developing  adaptive  methods  to  achieve  probabilistic  convergence  results  and  enhance  performance. 

The  purpose  of  this  STTR  Phase  I  is  to  demonstrate  the  feasibility  of  the  proposed  approach  through  its 
application  to  a  baseline  problem  involving  at  least  two  interactive  disciplines  and  two  levels  of  fidelity. 
Therefore,  objective  (5)  is  the  application  of  the  proposed  methods  to  a  flying  wing-type  UAV  design 
function  integrating  information  from  structural  and  fluid  models. 

4.  METHODS 

This  section  describes  in  detail  the  methods,  assumptions,  and  procedures  used  in  this  work.  We  begin 
(Section  4.1)  with  an  overview  of  the  optimization  process,  contrasting  the  traditional  approach  and  the 
multifidelity  framework  envisioned  in  this  work.  This  is  followed  (Section  4.2)  by  a  discussion  of 
surrogate  models,  including  the  rationale  for  our  choice  of  radial  basis  functions.  Section  4.3  describes 
the  optimization  methods,  and  Section  4.4  documents  the  computational  analyses  used  in  this  work. 

4.1  OVERVIEW 

In  the  conventional  view  of  the  design  process,  decisions  evolve  over  time.  At  the  conceptual  design 
stage  and,  increasingly,  in  preliminary  design,  a  large  array  of  options  and  cross-disciplinary  trades  are 
considered,  even  though  the  knowledge  of  the  system  is  necessarily  imprecise  (Lewis  and  Mistree,  1998). 
This  is  traditionally  the  design  space  exploration  phase,  where  individual  disciplines  must  be  able  to 
challenge  the  constraints  of  other  disciplines  (Holden  and  Keane,  2004)  in  order  to  come  up  with  effective 
solutions.  An  outcome  of  these  analyses  is,  traditionally,  to  narrow  down  the  design  space  to  smaller, 
more  manageable  portions  which  are  subsequently  investigated  in  more  detail.  The  intermediate  level  of 
fidelity,  though  more  expensive,  is  then  used  to  refine  the  analysis,  add  geometric  detail,  and  increase  the 
physical  fidelity  of  the  models  used,  thus  reducing  the  uncertainty  in  the  design.  Cross-disciplinary 
interactions  are  retained  in  the  MDO  methods  used,  but  a  smaller  portion  of  the  design  space  is 
considered,  to  make  this  tractable  in  spite  of  the  higher  costs.  Finally,  the  detailed  design  stage  is  used  to 
further  refine  the  analysis,  narrowing  the  design  to  a  handful  of  options. 


Figure  1.  Integrated  Hierarchical  Framework. 


In  contrast  to  the  conventional  view,  the  proposed  multifidelity  multidisciplinary  framework  (illustrated 
in  Figure  1)  is  a  true  integrated  hierarchical  model  which  takes  into  account  (merges)  the  information 
from  multiple  fidelity  levels  and  manages  this  information  in  order  to  achieve  optimal  sampling  of  the 
design  space,  both  in  terms  of  design  variables  and  level  of  fidelity.  The  right-hand  side  in  Figure  1 
represents  analysis  codes  spanning  three  hypothetical  disciplines,  e.g.,  fluids,  structures,  and  guidance  and 
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control.  Each  analysis  code  is  represented  by  a  O  symbol  and  is  ranked  according  to  its  level  of  fidelity, 
i.e.,  the  level  of  the  physics  included  in  the  model  and  the  level  of  geometric  detail.  Within  the 
multidisciplinary  multifidelity  framework,  codes  with  compatible  fidelity  levels  are  first  coupled  together 
using  cross-disciplinary  links  (CDL).  Such  groupings  are  the  staple  of  multidisciplinary  optimization 
(MDO)  and  include  not  only  the  code  associations  but  also  the  nature  of  the  coupling  between  them. 
Within  each  group,  not  only  are  variables  appropriately  exchanged  between  disciplines,  but  the  processes 
for  information  exchange/update  are  well-defined.  For  example,  the  disciplines  within  a  group  can  be 
tightly  integrated  or  loosely  coupled,  physical  interactions  accounted  for  explicitly  or  through 
subiterations  between  codes,  and  so  on.  The  left-hand  side  of  the  figure  symbolizes  the  two  key 
elements  which  underlie  the  proposed  approach.  The  first  is  the  merging  of  all  multifidelity  data  into  a 
global  probabilistic  response  surface-based  metamodel  (see  Section  4.2.2).  This  “Integrated  Hierarchical 
Metamodel”  builds  on  the  foundation  provided  by  multifidelity  radial  basis  function  (RBF)  data  fusion 
methods  (Reisenthel  et  al.,  2006a),  and  optimization  methods  using  probabilistic  RBF  networks  (see,  e.g., 
Reisenthel  and  Lesieutre,  2007).  The  second  element  and  a  key  addition  to  these  ideas  is  the  use  of  a 
Multifidelity  Sequential  Kriging  Optimization  (MFSKO)-like  “Expected  Improvement  Function” 

(Section  4.3.4)  to  manage  when  high-,  medium-,  and  low-fidelity  analysis  tools  are  appropriate  as  the 
“next  move”  in  design  space  exploration/exploitation  and  metamodel  updates.  The  goal  of  balancing 
exploration  and  exploitation  in  a  mathematically  rigorous  way  is  achieved  by  using  an  integrated  criterion 
which  determines  both  the  location  and  fidelity  level  of  subsequent  searches.  This  criterion  is  based  on 
Kennedy  and  O'Hagan's  (2000)  original  work,  and  subsequent  refinements  by  OSU  members  of  our  Team 
(Huang  et  al.,  2006,  Schenk  et  al.,  2005). 

The  work  presented  in  this  report  focuses  primarily  on  the  multifidelity  aspects  of  aerospace  vehicle 
analysis  and  design.  Thus,  the  “isofidelity”-groupings  depicted  in  Figure  1  represent  an  encapsulation  of 
existing  analysis  tools,  solution  methods,  and  information  exchange  processes  which  can  operate 
independently  of  each  other.  The  central  question  addressed  in  this  work  is  how  to  merge  and  manage  the 
information  from  these  multiple  fidelity  groupings. 

4.2  SURROGATE  MODELS 

At  the  heart  of  the  multifidelity  multidisciplinary  framework  depicted  in  Figure  1  reside  (1)  a 
dynamic  metamodel  and  (2)  an  intelligent  management  tool.  As  shown  in  Section  4.2.1,  the  choice  of 
surrogate  for  the  metamodel  has  a  number  of  consequences  which  must  be  taken  into  account. 

Section  4.2.2  details  the  multifidelity  radial  basis  function  surrogate  models  used  in  this  work. 

4.2.1  BACKGROUND/RATIONALE 

A  summary  comparing  the  pros  and  cons  of  four  different  classes  of  surrogate-based  methods  is  provided 
in  Table  1.  This  summary  organizes  information  gleaned  in  large  part  from  the  review  papers  of 
Forrester  et  al.  (2009)  and,  to  a  lesser  extent,  Queipo  et  al.  (2005),  and  supplemented  with  our  own 
experience.  The  table  is  organized  as  follows.  The  four  classes  of  methods  (polynomial,  radial  basis, 
Kriging,  and  support  vector  regression)  are  listed  across  the  top,  with  their  attributes  listed  by  row.  These 
key  attributes  are  categorized  according  to  (a)  the  forming  of  the  metamodel  surrogate,  (b)  considerations 
pertinent  to  optimization,  and  (c)  multifidelity  modeling.  While  many  of  the  table  entries  are  qualitative 
(e.g.,  “Low”  vs.  “Med.  Low”),  their  intent  is  to  give  a  sense  of  the  relative  performance  of  one  method  over 
another.  Whenever  possible,  this  sense  of  performance  is  quantified  in  terms  of  scaling  with  respect  to 
number  of  points,  search  space  dimensionality,  and  so  on.  Particular  method  variants  may  cross 
boundaries,  and  not  all  methods  are  represented,  so  this  is  not  intended  to  be  a  comprehensive  guide  but, 
rather,  a  reference  to  help  explain  where  the  proposed  methods  (highlighted)  fit  within  the  spectrum  of 
existing  approaches. 
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Table  1.  Characteristics  of  various  surrogate-based  methods  &  assessment  of  suitability  for  multifidelity  multidisciplinary  optimization. 


PRG 

MLS 

RBF 

qParam. 

RBF 

Param. 

RBF 

KRG 

Blind  KRG 

SVR 

Param. 

SVR 

Assumption  requirements 

High 

Medium 

Med.  Low 

Low 

Low 

V.  Low 

V.  Low 

Med.  Low 

V.  Low 

Landscape  complexity 

Low 

Medium 

Med.  High 

High 

High 

V.  High 

V.  High 

Med.  High 

V.  High 

Ability  to  incorporate  derivatives 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Unknown 

Unknown 

Generalization  error 

Medium 

Low 

Med.  Low 

Low 

Low 

V.  Low 

V.  Low 

Low 

V.  Low 

o 

Time  spent  training  or  retraining  the  metamodel 

Low 

High 

Low 

Low 

Medium 

High 

V.  High 

V.  High 

V.  High 

J2 

CO 

c 

O 

o 

a3 

Solution  effort  in  terms  of  number  of  parameters  searched 
(independently  of  GCV)* 

Number  of  potential  regressors,  given  m  data  points! 

0 

(m+d)l 
in  I  d ! 

0 

(m+d)l 

mid! 

0 

m 

0 

m 

1 

m 

2k 

m 

2/c2  +  1 

m 

2n 

m 

2  n  +  1 

m 

o 

5 

Maximum  number  of  data  points  m 

High 

V.  High 

<  1000 

<  1000 

<  1000 

<  500 

<  500 

High 

High 

-2 

CD 

Estimate#  of  maximum  design  space  dimensionality  k 

k<  20 

k<  20 

High 

High 

High 

k<  20 

k<  20 

High 

High 

Ease  of  implementation 

V.  High 

High 

High 

High 

Medium 

Medium 

Low 

Low 

Low 

Possibility  of  error  prediction  estimates 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Current  availability  of  error  estimates 

Unknown 

Unknown 

Yes 

Yes 

Unknown 

Yes 

Yes 

Yes 

Yes 

Adaptivity  (ability  to  fit  complex  local  behavior) 

Low 

Medium 

Medium 

High 

High 

High 

High 

Low 

High 

Cost  of  search 

Low 

High 

Low 

Low 

Low 

Low 

Low 

Low 

Low 

c: 

o 

CD 

N 

Basic  effort  associated  with  a  2-stage  infill,  in  terms  of 
number  of  parameters  searched  (independently  of  GCV)* 

0 

0 

0 

0 

1 

2k 

2  k2  +  1 

2  n 

2  n  +  1 

s 

-iS 

O' 

1 -stage11  infill  effort  in  terms  of  number  of  parameters 
searched  (independently  of  GCV)* 

Unknown 

Unknown 

k  +  1 

k  +  1 

k  +  2 

3  k+  1 

2k2  +  k  +  1 

2n  +  k  +  1 

2n  +  k  +  2 

Overall  suitability  for  optimization 

Med.  Low 

Low 

Medium 

High 

High 

Med.  High 

Medium 

Low 

Low 

Suitability  for  global  optimization 

V.  Low 

Low 

Medium 

High 

High 

Med.  High 

Medium 

Unknown 

Unknown 

Multifidelity 

Approximation  model 

or  surrogate  management  framework  used 

Affine  correction 
process 

Trust  region  methods 

• low  level  of  complexity 
between  fidelity  levels 
‘lack  of  flexibility  in 
terms  of  point/fidelity 
selection 

•Kennedy-  "Dimensionality 

O'Hagan  augmentation 

+  Least  Squares  Estimation 
•models  complex  relationships 
•flexible  point/fidelity  selection 
•fast 

Co-Kriging  (Kennedy- 
O’Hagan  +  Kriging) 
•potentially  high  level 
of  complexity 
between  fidelity  levels 
•flexible  point/fidelity 
selection 

Unknown 

Notes:  *  assumes  k-dimensional  design  space  and  n  support  vectors  #  Maximum  dimensionality  for  PRG  and  MLS  limited  due  to  large  number  of  regressors; 

t  assumes  m  data  points  and  polynomial  of  degree  d  maximum  dimensionality  for  KRG  and  Blind  KRG  limited  due  to  computational  effort; 

H  (conditional  lower  bound  approach);  estimated  scaling  potential  maximum  dimensionality  for  RJBF  variants  is  uncertain,  although  our  Team  has 

(except  KRG,  from  Forrester  et  al.,  2009).  successfully  tested  interpolating  RBF  solutions  of  random  data  in  100+  dimensions. 
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It  is  important  to  keep  in  mind  that  the  applicability  or  optimality  of  a  given  method  is  problem- 
dependent.  Thus,  for  example,  a  method  which  seems  unattractive  due  to  computational  demands  (e.g.,  if 
it  involves  layer  of  optimization  upon  layer  of  optimization)  may  be  worth  the  effort,  depending  on  how 
costly  a  single  high-fidelity  evaluation  is.  Paradoxically,  Forrester  et  al.  (2009)  note  that,  for  those  very 
applications,  i.e.,  those  most  likely  to  be  solved  in  a  highly  parallel  computing  environment,  “setting  up 
and  searching  a  surrogate  of  any  kind  can  (become  the)  bottleneck.  (...)  This  fact  will  always  limit  the 
amount  of  time  we  can  dedicate  to  the  building  and  searching  of  surrogates.” 

In  part  for  this  reason,  the  multifidelity  multidisciplinary  optimization  approach  recommended  by  our 
Team  centers  on  radial  basis  function  (RBF)  surrogation,  including  several  of  its  variants  (fixed  basis 
RBF,  quasi-parametric  RBF,  and  parametric  RBF).  As  depicted  in  Table  1,  this  class  of  methods  appears 
to  offer  qualities  which  very  nearly  approach  those  of  Kriging  (in  terms  of  modeling  ability,  flexibility, 
and  generalization  properties)  but  with  significant  advantages  in  terms  of  performance.  The  RBF  class  of 
surrogates  covers  a  wide  range  of  methods,  from  simple  fixed  basis  RBF  to  fully  parametric  RBF 
approaching  the  complexity  of  Kriging.  Within  this  class,  we  have  the  freedom  to  trade  generality  for 
performance,  depending  on  which  parameters  are  solved  by  optimization  (parametric  RBF)  versus  which 
parameters  are  fixed  (RBF).  It  is  argued  that  a  particularly  efficient  middle  ground  is  the  quasi- 
parametric  RBF  framework,  where  metaparameters  evolve  according  to  a  prescribed  algorithm.  Such 
algorithms  may  be  rooted  in  heuristics  and  are  designed  to  ensure  adaptivity  while  maintaining  high 
performance. 

4.2.2  MULTIFIDELITY  RADIAL  BASIS  FUNCTION  SURROGATE  MODELS 

In  this  section,  multifidelity  modeling  issues  are  described  for  models  involving  radial  basis  functions. 

We  describe  two  modeling  formulations  and  variance  and  bias  errors.  The  first  derives  from  previous 
research  in  which  radial  basis  functions  model  the  systematic  errors  of  each  surrogate  system  in  relation 
to  higher  fidelity  systems.  The  second  formulation  treats  fidelity  as  a  design  dimension  and  adjusts  the 
specific  scaling  in  the  radial  function. 

4.2.2.1  Modified  Kennedy  and  O’Hagan  Scheme 

The  first  and  primary  modeling  formulation  is  adapted  from  Kennedy  and  O’Flagan  (2000)  involving  both 
linear  and  Kriging  models.  This  was  the  approach  used  in  Huang,  Allen,  Notz,  and  Miller  (2006).  The 
“cascading”  formulation  is: 

f(x)=Mx)  +  d,(x)  (/  =  2, 3,  m)  (1) 

where  fi  is  the  function  evaluated  at  fidelity  level  /  and  r)/(x)  is  independent  of  /i(x),  /2(x),  .  ,.,f. i(x)  and  x 
is  a  d  dimensional  decision  vector.  Hypothetically,  different  levels  of  /  also  could  be  different  levels  of  a 
discrete  factor  for  the  same  fidelity  level.  For  convenience  in  the  notations  below,  we  also  let: 

/i(x)  =  <5i(x).  (2) 

Note  that  for  /  =  2,  3,  . . .,  m,  d /(x)  can  be  understood  as  the  “systematic  error”  of  a  lower-fidelity  system, 
(/— 1 ),  as  compared  to  the  next  higher- fidelity  system,  /.  In  these  cases,  <5/(x)  is  usually  small  in  scale  as 
compared  to  f(x),  otherwise  there  will  be  no  reason  for  the  lower-fidelity  system  to  exist.  Of  course,  in 
some  cases,  physics  may  be  missed  in  low-fidelity  model(s)  and  d/(x)  will  be  large. 

In  radial  basis  metamodeling,  the  response  is  assumed  to  be  a  linear  model.  We  use  another  radial  basis 
function  to  model  the  departure  of  the  lowest-fidelity  system,  A(x),  as  well  as  the  difference  between 
systems,  A(x)  (/  =  2,  3,  . . .,  m).  Therefore,  we  have 

S{x,Si, . . . ,Sm)  =  b(\,S i , . . . ,<S))TP/  +  £/  (/=  1,  2,  m)  (3) 

where  b/  and  P/  are  the  basis  functions  and  coefficients,  respectively,  of  a  linear  model.  Also,  the  sets 
Si,..  ,,Sm  contain  the  input  points  for  selected  past  runs  at  the  m  levels  of  fidelity.  In  general,  not  all 
previous  runs  are  included  in  the  model  so  that  there  can  be  lack  of  fit  estimation. 
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The  linear  model  portion  is  richer  in  Eq.  (3)  than  in  Huang,  Allen,  Notz,  and  Miller  (2006)  because  it 
contains  dependence  on  design  points.  The  basis  functions  of  the  linear  model  here  are  radial  centered  on 
the  design  points.  By  basing  the  model  centering  on  design  points  there  is  no  nonlinearity  associated  with 
picking  runs. 

In  our  implementations,  we  use  Gaussian  basis  functions  or,  alternatively,  reciprocal  multiquadrics 
centering  in  both  cases  on  k  selected  runs.  In  this  report,  we  focus  somewhat  arbitrarily  on  Gaussian  basis 
functions  but  our  experiences  with  the  use  of  multiquadrics  were  similar.  For  a  single  fidelity  level  the 
Gaussian  basis  function  is: 

cHx,x/)  =  exp[-7(x  -  xj)'(x  -  x,)]  (j  =1,2,...,  k).  (4) 

where  y  is  an  adjustable  scale  parameter.  In  this  report,  we  focus  on  y  as  a  fixed  parameter.  Yet,  the 
adjustment  of  y  as  the  algorithm  evolves  is  considered  promising  for  further  improvements  to  numerical 
results  on  test  problems.  In  general,  our  Team  feels  that  desirable  algorithms  that  increase  y  as  the 
algorithm  are  most  likely  to  foster  improved  results. 

Therefore,  the  general  form  of  the  prediction  model  at  the  point  x  at  fidelity  level  /  is: 

fl(x)  =  =  S/CS!  <Kxwx;)A,/  .  (5) 


With  these  assumptions,  the  linear  model  design  matrix,  X,  is  somewhat  sparse.  For  example,  with  three 
fidelity  levels  and  k\,  k2,  and  ki  runs  included  in  the  model  and  nl,  m,  and  773  total  runs  at  each  of  three 
fidelity  levels,  the  X  matrix  is: 


4>(Xl,X1) 

-  <Kx1(xki) 

0 

0 

^(Xn^Xi) 

-  ^(Xn^X^) 

4>(x„1+i,x1) 

-  4>(x„1+i,xki) 

4>(x„1+l,Xkl+1) 

4>(xn1+l»Xkl+kl) 

0 

‘KXnj+ns.M) 

-  ‘Kx^+n^X^) 

^(.Xnj+m.X^+i) 

■"  ‘KXni+nj/Xlij+lt,) 

<t>(Xn1+nI+l.X1) 

4>(xn1+nI+l<Xkl) 

<Kxr!l+nI+l>Xkl+i) 

<KxrIl+n2+l<Xkl+kl) 

... 

-4KXni+n2+nslXi) 

4*  xki) 

0  (xni+ni+n3 ,  xki+i) 

4^(xr,1+ni4-n3,Xkl4-k2) 

— 

In  our  numerical  results,  we  fitted  Eq.  (5)  by  solving  the  normal  equations  using  singular  value 
decomposition  (SVD)  and  ANSI  C  code.  This  code  was  significantly  modified  from  code  in  Press, 
Teukolsky,  Vetterling,  and  Flannery  (2007).  The  modifications  include  facilitating  multivariate 
modeling,  covariance  estimation,  and  error  prediction. 

The  modified  version  is  order  (ni  +. . .+  nmf  but  it  requires  generally  negligible  computing  time  compared 
with  likelihood  optimization  in  sequential  Kriging  optimization  (Huang,  Allen,  Notz,  and  Miller  2006). 
Also,  linear  model  estimation  using  SVD  is  generally  more  reliable  and  reproducible  in  part  because  of 
limited  sensitivity  to  numerical  issues. 

4.2.2.2  An  Alternative  Augmented  Dimensionality  Formulation 

Our  team  also  considered  alternative  modeling  formulations  for  addressing  fidelity  based  on  scaling 
inputs  (Reisenthel  et  al.,  2006).  In  some  of  these  formulations,  global  adjustments  to  distances  are  made 
based  on  the  differences  between  fidelity  levels  of  the  runs.  Assume  that  f(z)  is  a  1  x  m  vector  with  a  1 
when  i  is  the  fidelity  level  of  run  i  or  0  otherwise.  In  this  approach,  we  consider  adjustable  parameters, 
y, >  0  ,  for  each  of  the  design  dimensions  i=\,...,d  and  Gij>  0  for  all  fidelity  levels  j  =  1  ,...m  and 
i  =  1,. .  .m.  The  revised  density  is: 

^(Xl,f(0,X;,f(/-))  =  exp{-Xfc=1  F,(xk.i  -  xk]y.  G,||f(0  -  f(/)  II2}  (6) 
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where  x*,,-  is  the  klh  element  of  the  vector  x„  The  model  is  simply: 

fix,q)  =  E/= l . » <t>[x,<7,x7,f0)]/?7 


And  the  design  matrix  is: 

<t>(xn»rn»xl<rl) 


<J)(xl»rl»xn>rn) 

d>(x„,rn,xn,rn). 


(7) 


where  we  have  assumed  that  all  runs  have  been  included  in  the  model  as  an  example.  Therefore,  Eq.  (6) 
would  represent  a  saturated  model  that  generally  passes  through  all  of  the  design  points. 

This  alternative,  more  concise  formulation  has  the  disadvantage  that  the  Gy  weight  parameters  must  be 
adjusted  heuristically.  Also,  there  is  no  guarantee  that  the  process  will  converge  to  accurate  models  of  the 
bias  or  systematic  errors  or  the  predictions  of  the  highest  fidelity  model  will  converge.  Yet,  the 
augmented  dimensionality  approach  includes  far  fewer  parameters  to  be  estimated  and  can  be  expected  to 
offer  efficiency  advantages  for  cases  in  which  assumption  parameters  match  problem  properties  closely. 

4.2.2.3  Modeling  Prediction  Variance 

A  possible  advantage  of  Kriging  models  over  radial  basis  functions  is  that  they  provide  a  model  of  mean 
and  uncertainty  even  if  the  model  passes  through  all  the  design  points.  Kennedy  and  O’Hagan  (2000)  and 
others  have  clarified  the  Bayesian  nature  of  the  intervals.  Yet,  the  interpretation  and  validity  of  the 
derived  intervals  in  realistic  cases  is  not  fully  understood.  At  the  same  time,  the  computational 
performance  in  test  problems  of  methods  based  on  these  Kriging  error  estimates  appears  to  be  good. 

By  contrast,  the  radial  basis  functions  considered  in  previous  sections  are  associated  with  regression  type 
intervals.  If  the  prediction  model  passes  through  all  the  design  points,  the  standard  regression  intervals 
indicate  zero  prediction  errors  at  all  points.  This  follows  because  the  standard  intervals  are  based  on  so- 
called  variance  errors  and  therefore  derive  entirely  from  the  assumed  repeatability  or  random  errors.  We 
consider  four  possible  approaches  to  address  predictions  of  prediction  uncertainty: 

1 .  Assumed  random  errors  -  Some  of  the  computational  results  from  the  team  derived  from  ad  hoc 
assumptions  of  random  error  standard  deviation,  a,  inserted  into  the  prediction  variance  formula: 

o2l(x)=  Var  \y prediction^, l)\  =  c02b(x,su.  ■  .,Si)T(X'X)~  lb(x,Si,...,Si).  (8) 

This  approach  has  some  intellectual  coherence  in  the  context  of  FEM.  This  follows  because  there 
is  often  a  known  or  typical  uncertainty  associated  with  computer  code,  e.g.,  related  to  the  mesh 
size.  Yet,  the  errors  are  not  like  typical  repeatability  errors. 

2.  Simply  estimated  random  errors  -  A  selected  set  of  runs  can  be  dropped  from  the  model  and  used 
for  estimating  o  as  is  common  in  regression.  In  our  computational  results,  we  simply  drop  the 
most  recently  collected  r  runs  from  fitting.  The  standard  Analysis  of  Variance  estimate  of  the 
random  error  is  then: 

o-o  =  V[(I"- Xpf][(Y- Xp)]/(n-r)  (9) 

which  is  based  on  an  equal  variance  assumption  as  is  standard  in  regression. 

3.  PRESS  estimated  random  errors  -  The  standard  regression  cross  validation  method  of  leaving  out 
each  observation,  fitting  with  the  others,  and  then  rotating  can  be  applied  to  estimate  c.  This 
estimate  is  based  on  an  average  of  calculations  from  Eq.  (9). 

4.  Bias  plus  variance  -  Research  by  the  OSU  members  of  our  Team  has  recently  invented  formulas 
for  the  variance  and  the  bias.  Such  approaches  are  well-founded  and  make  reference  to  high- 
order  polynomial  true  models.  They  offer  an  intellectually  coherent  way  to  assign  uncertainty 
even  for  cases  like  FEM  with  known  perfect  repeatability.  Tseng  (2007)  provides  several  useful 
formulas. 
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Option  4  is  the  most  attractive.  However,  in  practice  estimating  bias  errors  requires  the  assumption  of  the 
model  form  of  the  true  model,  e.g.,  a  fourth-order  polynomial,  and  a  distribution  of  the  unknown 
coefficients.  The  associated  arbitrariness  of  related  assumptions  makes  these  methods  seem  less 
attractive.  Yet,  Kriging  models  have  associated  ambiguities  also  relating  to  the  validity  and  convergence 
of  the  assumed  correlation  parameters.  At  present,  we  focus  on  options  (1)  and  (2)  from  the  above  list 
and  leave  options  (3)  and  (4)  for  future  research. 

4.3  OPTIMIZATION  METHODS 

This  section  describes  the  central  optimization  problem  and  solution  method  framework.  The 
framework  is  based  on  sequential  optimization  using  radial  basis  function  metamodels.  It  is  based  on 
a  search  after  every  computational  analysis4  for  the  location  of  the  next  run  in  the  parameter  and 
fidelity  space  that  maximizes  the  expected  improvement  in  a  manner  similar  to  Huang,  Allen,  Notz, 
and  Miller  (2006).  Next,  the  expected  improvement  function  is  described  and  method  variations  based 
on  it  are  defined.  Finally,  the  developed  Visual  Basic,  ANSI  C,  and  Octave  codes  developed  in  this 
work  are  briefly  characterized. 

4.3.1  THE  OPTIMIZATION  PROBLEM 

Suppose  there  are  a  total  of  m  systems  to  draw  evaluations  from,  including  the  real  and  the  surrogates. 
Denote  the  output  functions  of  these  systems  in  increasing  order  of  fidelity  by  /I(x),/2(x),  ...,/„(x), 
where  x  is  the  input  vector.  Therefore,  f(\)  has  the  lowest  fidelity,  and  fjx)  has  the  highest  fidelity. 

As  mentioned  previously,  the  highest-fidelity  system  is  the  system  of  interest,  therefore  the  goal  is  to 
minimize  /,„(x)  within  the  feasible  region,  y,  i.e. 

min /,»(*)  .  (10) 

xe/  v  7 

We  consider  the  systems  as  black  boxes  that  provide  no  information  other  than  measurements  of  the 

outputs.  Denoting  by  d  the  dimension  of  the  input  space,  we  assume  that  the  feasible  region  X  c  Rd  is 

connected  and  compact. 

Each  system  is  associated  with  a  cost-per-evaluation,  which  is  denoted  by  Cy,  C2,  ...,  C„„  respectively.  In 
this  research,  the  total  cost  of  all  evaluations  measures  the  efficiency  of  the  optimization  scheme. 

Usually,  a  lower-fidelity  evaluation  is  cheaper  than  a  higher-fidelity  evaluation,  i.e.,  Cy  <  C_?  <  ...  <  C,„. 
Also,  for  now  we  assume  that  the  cost  of  even  the  cheapest  system  is  somewhat  expensive,  such  that  it  is 
“worthwhile”  to  regenerate  a  radial  basis  function  metamodel  in  order  to  determine  the  next  search 
location. 

In  addition,  the  measurements  of  a  system  output  may  contain  random  error  or  noise.  For  each  system, 
we  assume  that  random  errors  from  successive  measurements  are  independent  identically  distributed 
(IID)  normal  deviates. 

4.3.2  GENERAL  OPTIMIZATION  FRAMEWORK 

The  outline  for  the  proposed  Multiple  Fidelity  Sequential  Radial  Basis  Optimization  (MFSRBO)  is  as 
follows: 

Step  1  :  Initialize  radial  basis  function  parameter(s),  e.g.,  y  and  the  residual  parameters,  e.g.,  r,  the 
number  of  computational  analyses  for  residual  or  error  estimation. 

Step  2  :  Perform  an  initial  experimental  design  involving  all  fidelity  levels.  Section  4.3.3  describes  the 
options  and  focuses  on  n  run  Latin  hypercube  initial  experiments  following  Huang,  Allen,  Notz, 
and  Miller  (2006). 


4  also  referred  herein  as  “experimental  run” 
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Step  3  :  Fit  the  radial  basis  function  of  the  selected  type  (see  Section  4.2.2)  to  all  available  data  using 
singular  value  decomposition  (SVD)-based  least  squares  estimation. 

Step  4:  Find  the  location  and  fidelity  level  of  the  new  evaluation  that  maximize  the  augmented  expected 
improvement  (El)  function.  If  the  El  function  is  sufficiently  small,  go  to  Step  6. 

Step  5:  Conduct  an  evaluation  at  the  selected  point  from  Step  4.  Go  to  Step  3. 

Step  6:  Perform  an  additional  search  and/or  evaluation  to  evaluate  solution  quality.  For  example,  apply 
trust  region  evaluation  which  includes  a  check  for  Karush-Kuhn-Tucker  point  convergence. 

Note  the  removal  of  a  diagnostic  step  from  the  Fluang,  Allen,  Notz,  and  Miller  (2006)  because  there  is  no 
need  to  assume  zero  expected  mean  for  the  biasing  functions,  8/.  This  constitutes  one  advantage  of 
sequential  radial  basis  optimization  (MFSRBO)  methods  proposed  here  and  may  be  expected  to  lead  to 
relatively  robust  performance,  particularly  in  cases  involving  large  systematic  errors. 

By  convention,  Step  1  is  also  referred  to  as  the  “initial  fit”  stage,  while  Steps  4  and  5  are  called  the 
“infill”  or  “update”  stage.  The  sequentially  added  evaluations  are  also  called  the  “infill”  or  “update” 
points.  Note  that  the  proposed  method  differs  from  its  predecessors  mainly  in  the  following  two  aspects: 
1 )  the  radial  basis  function  metamodel  is  generated  using  multiple  fidelity  data,  and  2)  the  El  formulation 
takes  into  account  not  only  the  location  but  also  the  fidelity  level  of  an  infill  point. 

4.3.3  DESIGN  FOR  INITIAL  FIT 

Step  2  in  the  previous  section  involves  an  experimental  design  for  an  initial  fit.  This  is  essentially  the 
same  as  the  initial  step  in  methods  based  on  Kriging  models  (e.g.,  Fluang,  Allen,  Notz,  and  Miller 
2006).  Again,  assume  the  number  of  factors  is  d.  The  methods  in  Huang,  Allen,  Notz,  and  Miller 
used  an  n  =  10  x  d  run  maximum  minimum  distance  Latin  hypercube  (from  the  Ihsdesign  function  in 
MATLAB®)  for  the  lowest  fidelity  system  and  then  selected  a  subset  for  higher  fidelity  systems.  If 
the  system  includes  random  errors,  d  replicates  are  added  to  all  levels  at  locations  where  the  best 
responses  are  found. 

As  another  alternative  for  future  exploration,  we  suggest  an  initial  design  that  includes  all  points  on  a 
grid  for  the  lowest  fidelity  level.  For  each  higher  fidelity  system,  a  subset  of  the  previous  points  is 
selected  including  the  q%  of  the  best  results  from  lower  fidelity  levels.  In  most  of  the  numerical 
examples  shown  in  this  report,  q  =  25%.  The  rational  for  using  a  grid  instead  of  a  Latin  hypercube  is 
that  the  projection  properties  of  Latin  hypercubes  are  likely  not  relevant  since  all  systems  are 
generally  important.  Also,  radial  basis  functions  are  not  affected  by  the  numerical  errors  for  Kriging 
models  associated  with  equally  spaced  inputs. 


4.3.4  EXPECTED  IMPROVEMENT  FUNCTIONS 


The  expected  improvement  function  fits  into  the  search  framework  in  Section  4.3.2  and  is  used  to 
select  the  next  “infill”  point.  The  function  balances  the  desire  to  improve  the  search  criterion  with  the 
desire  to  reduce  uncertainty  in  a  manner  reminiscent  of  trust  region  point  selection.  The  methods 
explored  here  are  a  straightforward  application  of  the  formulation  in  Huang,  Allen,  Notz,  and  Miller 
(2006)  except  that  they  are  based  on  radial  basis  functions  instead  of  Kriging  models. 


Therefore,  for  Multiple  Fidelity  Sequential  Radial  Basis  Optimization  (MFSRBO),  we  propose  the 
following  augmented  Expected  Improvement  function: 


EI{x,l)  =  E  max  fm(  x  )-fm(x)-am(x)e,0\  ax(x  ,l)-a2(x,  l)-a3(l) 


(ID 


where  £  is  normally  distributed  N[0,1]  and 

al(x,l)  =  function  that  discounts  systems  relating  to  their  predicted  accuracy,  (12) 
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a 


and 


,(x,l)=  (l  -  a0  / ^Jdf(x)+b jj) 
Cost  , 


(13) 


Of,  (  /  ) 


Cost, 


If  the  problem  is  deterministic,  i.e.,  with  zero  repeatability  error,  then  we  (generally)  use  (XiixJ)  =  1  and 
J  m\X  )  is  the  minimum  of  the  already  sampled  points.  More  generally,  if  at  least  a  single  system  has 
noise,  in  Eq.  (11),  x*  stands  for  the  current  “effective  best  solution”  defined  by 


* 

x  = 


argmax 

xe[x1,x2,...xn\ 


ll(  x) 


(14) 


where  u  ( x )  =  -  fm{  x )  -  p  bj  x )  .  The  formula  from  Jones,  Schonlau,  and  Welch  ( 1 998)  is: 
m^x[f(x*)-f(x)-o{x)£,Q\ 


- 

f 

J  m 

v  )~fjx) 

o(x ) 

\ 

m  v  /  / 

+  °m(x)<l> 


fm(X 


°m\X> 


(15) 


where  ip  is  the  cumulative  normal  and  </>  is  the  normal  density.  We  need  the  two  formulas  for  the 
variance,  i.e.,  with  and  without  bias.  In  Huang,  Allen,  Notz,  and  Miller  (2006),  oq  was  chosen  as 


a 


1(x,/)=corr  fm(x)+£1am(x),fl(x)+e2cr/(x) 


(16) 


where  “corr“  is  the  correlation  function  and  ei  and  A  are  independent  normally  distributed  deviates. 

This  choice  had  the  following  desirable  properties: 

1 .  It  is  easy  to  compute  in  the  context  of  Kriging  models, 

2.  When  /  =  m,  it  equals  1 . 

3.  For  deterministic  systems,  it  equals  zero  at  past  design  points  since  G/(x)  is  zero. 

Unfortunately,  for  radial  basis  function  systems,  Eq.  (16)  is  not  easy  to  compute.  In  the  “augmented 
dimensionality”  version  of  our  method,  we  inserted  an  estimate  of  Eq.  (16)  based  on  numerical 
integration  (see  Appendix  A).  For  other  methods,  we  used  the  following  formulation  which  is  relatively 
easy  to  compute: 


oq  (x,  / )  = 


°i  (xj 


l/mU)-//(x)|  +  am(x) 


(17) 


where  ||  is  the  absolute  value.  This  approach  has  the  above-mentioned  desirable  properties  in  the 
context  of  radial  basis  functions. 


4.3.5  METHODS  VARIATIONS 

Every  combination  of  modeling  method,  expected  improvement  function  (assumed  random  error, 
simply  estimated  random  errors,  PRESS,  and  bias  plus  variance),  expected  improvement  function 
(deterministic  or  stochastic),  and  termination  sequence  (simple  or  trust  region)  hypothetically 
corresponds  to  a  method  variation.  Here,  we  consider  only  a  select  few  combinations,  characteristics 
of  which  are  summarized  in  Table  2. 

If  the  method  has  simple  termination  then  it  is  called  multifidelity  sequential  radial  basis  optimization 
(MFSRBO)  or,  for  cases  involving  only  a  single  level  of  fidelity,  sequential  radial  basis  optimization 
(SRBO).  If  the  method  terminates  by  applying  a  downhill  search  such  as  trust  region  augmented 
Lagrangian  methods,  then  we  refer  to  the  method  as  a  “hybrid”. 
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Table  2.  Variants  considered  in  the  computational  and  theoretical  approach. 


Method  Variant 

Prediction  Uncertainty 

El  Function 

Termination 

1.  MFSRBO:  augmented 
dimensionality 

Variance  Intervals  based  on 
assumed  prediction  errors 

Deterministic 

Simple 

2.  MFSRBO:  cascading 
model 

Variance  intervals  from  estimated 
random  errors 

Deterministic 

Simple 

3.  SRBO:  cascading  model 

Variance  intervals  from  estimated 
random  errors 

Stochastic 

Simple 

4.  Hybrid 

Variance  intervals  from  estimated 
random  errors 

Both 

Trust  Region  Augmented 
Lagrangian 

4.3.6  COMPUTER  CODES 

The  method  “MFSRBO:  augmented  dimensionality”  was  used  to  demonstrate  multifidelity  radial 
basis  function  optimization  on  a  real-world  aerospace  example  (Section  5.2.3).  This  method  was 
implemented  using  NEAR  RS,  a  radial  basis  function-based  regression  tool  whose  response  surface 
search  module  was  modified  to  accommodate  the  computation  of  the  multifidelity  expected 
improvement  function  as  the  optimization  driver.  In  NEAR  RS  the  multifidelity  cumulative  metamodel 
concept  relies  formally  on  a  dimensionality  augmentation  of  the  input  data  (control  variables)  space  by 
including  the  fidelity  level  as  an  additional  variable,  Eq.  (6).  Search  projections  are  required  to  counter 
this  dimensionality  augmentation  and  confine  the  optimization  to  the  subspace  corresponding  to  the 
desired  level  of  fidelity.  In  other  words,  we  seek  to  optimize  a  function  (or  its  metamodel 
representation)  on  the  high-fidelity  manifold,  regardless  of  where  low-fidelity  optima  may  be.  The 
point  of  view  adopted  here  is  that  the  lower  fidelity  data  act  as  “support”  vectors  aiding  the 
characterization  of  the  high-fidelity  space. 

NEAR  RS's  response  surface  search  module  is  written  in  GNU  Octave,  an  interpreter-based,  Matlab®- 
like  program  for  high  level  computation  that  is  freely  available  under  GNU  General  Public  License. 
While  very  versatile,  the  associated  code  was  not  judged  to  be  computationally  efficient  enough  for 
widespread  use.  This  is  due  in  part  to  the  nature  of  Octave,  as  opposed  to  the  execution  speed  of  a 
compiled  language,  such  as  ANSI  C,  and  in  part  to  the  cost  of  computing  numerically  the  correlation 
term  (Appendix  A)  at  each  evaluation  of  the  surrogate.  Additional  complications  include  the  use  of 
Options  1  and/or  2  (Section  4.2.2. 3)  for  the  modeling  prediction  variance,  ascribing  an  inputed 
variance  at  each  data  point  (“support”  vector)  according  to  what  is  believed  to  be  a  known  or  typical 
uncertainty  associated  with  a  particular  code  for  a  particular  problem.  This  approach  lacks  rigor, 
because  it  ascribes  an  essentially  epistemic  type  of  uncertainty  (i.e.,  no  assumptions  in  regard  to  the 
distribution  of  the  error  within  the  interval),  and  subsequently  treats  this  uncertainty  as  if  it  were 
aleatory  uncertainty5.  Furthermore,  what  could  be  principally  bias  error  is  treated  as  variance. 

Thus,  while  the  above  methods  have  enjoyed  empirical  success,  the  assumption  of  the  variance 
complicates  comparison  with  methods  that  are  not  based  on  assumed  random  errors  in  Huang,  Allen, 
Notz,  and  Zheng  (2005)  and  Huang,  Allen,  Notz,  and  Miller  (2006).  Also,  these  alternative  methods 
have  explicit  models  for  the  systematic  errors  through  their  “cascading”  model  formulation  and  thus 
likely  converge  to  a  true  picture  of  the  fidelity  hierarchy.  For  this  reason,  the  comparison  results 
presented  in  Section  5.2.2  are  based  on  method  variants  2-3  (SRBO:  cascading  models  and  MFSRBO: 
cascading  models)  and  were  implemented  in  a  new  code,  the  features  of  which  include: 

[1]  The  ability  to  call  Windows  simulation  and  read  results  using  shell  and  wait, 

[2]  Enumeration-based  (inefficient  but  repeatable)  expected  improvement  function  optimization 
(every  iteration,  the  code  evaluates  all  alternatives  on  a  grid  to  select  the  point  that  maximizes 
the  expected  improvement),  and 

[3]  Microsoft  excel  interface. 


5  In  other  words,  the  input  uncertainty  internal  is  treated  as  the  standard  error  of  a  stochastic  process. 
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The  code  developed  for  these  purposes  is  a  combination  of  ANSI  C++  and  Visual  Basic  (VB).  The 
developed  code  is  not  owned  or  copyrighted  by  any  participants  not  on  the  Team.  The  singular  value 
decomposition  (SVD)  estimation  code  was  based  on  numerical  recipes  code  which  was  significantly 
altered.  The  hybrid  designs  considered  here  are  only  for  the  convergence  results  established  in  the  next 
section.  Additional  hybrid  designs  are  proposed  for  future  research  to  offer  improved  computational 
results  while  achieving  rigorous  convergence. 

4.4  COMPUTATIONAL  ANALYSIS  METHODS 

The  computational  results  presented  in  this  report  (Section  5.2)  consist  of  several  test  problems  (e.g., 
Branin  function,  multifidelity  Rosenbrock  function)  and  two  applications  involving  computational 
structural  and  computational  fluid  modeling. 

For  the  structural  analysis  tool  we  used  McIntosh  Structural  Dynamics'  finite  element  code  CNEVAL 
(Lesieutre  et  al.,  1997).  The  wing  structure  was  modeled  using  conformal  quadrilateral  plate  elements 
cantilevered  at  the  root  chord.  Given  the  normal  aerodynamic  forces,  CNEVAL  provides  the  structural 
displacements,  stresses,  frequencies,  and  overall  structural  weight.  CNEVAL  was  used  to  predict  the 
static  aeroelastic  deformations  and  to  verify  that  structural  limits  were  not  exceeded.  CNEVAL  uses  the 
current  wing  planform  information,6  along  with  the  planfonn  thickness  information  and  material 
properties,  and  updates  the  finite-element  model  of  the  wing  which  is  subdivided  into  quadrilateral  panels 
used  to  distribute  the  aerodynamic  loads  for  the  wing  displacement  calculations.  Each  panel  is  divided 
into  2  triangles.  An  option  also  exists  to  further  subdivide  a  quad  element  into  4  triangles  for  greater 
accuracy.  Both  high-  and  low-fidelity  structural  modeling  were  carried  out  using  CNEVAL  by  changing 
the  finite  element  discretization,  as  described  in  Table  3. 


Table  3.  Finite  element  discretization. 


Wing  Design 

Fidelity 

Number  of  Structural  Elements 

Problem 

Level 

Chordwise 

Spanwise 

Unswept  tapered 

High 

10 

15 

4  triangle  elements  per  structural  guad 

Low 

2 

2 

2  triangle  elements  per  structural  guad 

UAV,  swept 

High 

10 

15 

2  triangle  elements  per  structural  guad 

Low 

4 

7 

2  triangle  elements  per  structural  guad 

A  version  of  NEAR's  intermediate-level  aerodynamic  prediction  code  M1SDL  (Lesieutre  et  al.,  1998, 
Dillenius  et  al,  1999)  was  used  as  the  high-fidelity  fluid  model.  MISDL  is  based  on  panel  methods, 
vortex  lattice  with  compressibility  correction  for  subsonic  flow,  and  other  classical  singularity  methods 
enhanced  with  models  for  nonlinear  vortical  effects.  The  method  is  applicable  to  subsonic  and  supersonic 
flight  vehicles  including  aircraft,  UAVs,  missiles,  and  rockets.  Circular  and  noncircular  cross  section 
bodies  are  modeled  by  either  subsonic  or  supersonic  sources/sinks  and  doublets  for  volume  and  angle  of 
attack  effects,  respectively.  Conformal  mapping  techniques  are  used  for  noncircular  bodies.  The  fm/wing 
sections  are  modeled  by  a  horseshoe-vortex  panel  method  for  subsonic  flow  and  by  first-order  constant 
pressure  panels  for  supersonic  flow.  Up  to  three  lifting  surface  sections  can  be  modeled  and  nonlinear  fin 
and  body  vortices  are  modeled.  MISDL  predicts  overall  aerodynamic  performance  coefficients  and 
detailed  aerodynamic  loading  distributions.  The  code  has  seen  extensive  use  in  missile  aerodynamic 
analysis  and  design,  including  aerodynamic  shape  optimization  and  multidisciplinary  design  optimization 
when  coupled  to  a  structural  finite  element  model.  In  recent  years,  MISDL  has  seen  use  in  the  analysis  of 
manned  and  unmanned  (UAV)  aircraft. 

OPTMIS  is  the  name  used  to  designate  the  version  of  MISDL  that  is  coupled  with  CNEVAL.7  The  cross- 
disciplinary  coupling  is  handled  by  the  OPTMIS  logic  through  successive  iteration  between  codes. 
OPTMIS  handles  interpolation,  displacement  and  load  transfer  operations  and  iterates  between 


6  i.e.,  the  (x,y)-projected  shape  of  the  wing 

7  OPTMIS  also  contains  an  optimization  module  of  its  own,  based  on  Powell's  conjugate  gradient  method,  but  the 
optimization  component  was  not  used  in  this  study. 
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aerodynamic  and  structural  predictions  until  a  prescribed  convergence  level  is  reached.  Convergence  is 
obtained  when  the  change  in  the  displacements  between  successive  iterations  is  less  than  a  user-prescribed 
threshold  at  all  structural  node  points.  For  the  purposes  of  this  study,  OPTMIS  was  enhanced  with  the 
computation  of  the  induced  drag.  Rather  than  incorporating  a  Trefftz  plane  analysis,  such  as  might  be 
done  in  an  Euler  code,  the  induced  drag  was  computed  using  a  direct  method,  as  follows. 

To  compute  the  induced  forces  on  the  wing,  the  Kutta-Joukowski  theorem  is  used:  F ,•=  pUjXF {  .  This 
induced  force  vector  is  computed  on  each  panel  as  the  cross  product  of  the  total  induced  velocity  on  the 
panel  (sum  of  all  vortex  lattice  induced  velocities;  the  freestream  velocity  being  omitted)  and  the  panel 
bound  vortex  vector.  The  drag  component  of  F .  from  each  panel  is  computed  and  summed  to  obtain  the 
total  induced  drag.  Several  test  cases  were  used  to  verify  the  calculation,  and  the  theoretical  limit  of 
Cd:  =  Cl2/(tiAR)  was  obtained  for  elliptically  loaded  wings.  An  option  also  exists  to  prescribe  a  cosine¬ 
shaped  spanwise  layout  of  the  panels8  to  increase  computational  accuracy.  The  low- fidelity  fluid 
modeling  amounted  to  simple  lifting  line  theory,  conveniently  implemented  within  OPTMIS  using  a 
degenerate  panel  layout,  as  shown  in  Table  4.  Code  robustness  issues  precluded  the  use  of  this 
implementation  for  the  case  of  the  UAV  wing  planform  and,  so,  for  this  case,  the  low-fidelity  model 
consisted  of  a  coarse  panel  discretization  aerodynamic  model,  the  details  of  which  are  also  furnished  in 
Table  4. 


Table  4.  Aerodynamic  Panel  Layout. 


Wing  Design 
Problem 

Fidelity 

Level 

Number  ofAerodynamic  Panels 
Chordwise  Spanwise 

Spanwise  Panel  Layout 

Unswept  tapered 

High 

10 

15 

Cosine  distribution 

Low 

1 

2 

Uniform 

UAV,  swept 

High 

10 

15 

Uniform 

Low 

4 

7 

Uniform 

The  isotropic  material  used  for  the  structure  was  aluminum.  We  used  the  following  material  constants: 


Table  5.  Material  Properties. 


Material  Property 

Value 

Young's  modulus 

10.5  X  106  lbf/in 2 

Poisson's  ratio 

0.3 

Density 

0.1  lb/in3 

Maximum  allowed  Von  Mises  stress 

22.5  X  103  lbf/in2 

The  structural  cross  section  geometry  was  self-similar  as  a  function  of  span,  and  idealized  as  shown  in 
Figure  2.  All  parameters  in  Figure  2  were  modeled  as  being  linearly  varying  from  root  to  tip,  with 
maximum  and  minimum  values  indicated  in  Table  6. 


Table  6.  Geometric  Parameters  (in  ft.). 


Figure  2.  Structural  cross-section. 


Geometric  Parameters 

Tapered  wing 

UAV  wing 

o 

o 

X 

0.2 

0.348 

X2_root 

0.3 

0 

TH_root 

variable 

0.025 

Xltip 

0.2 

0.103 

X2 _ tip 

0.3 

0 

TH_tip 

0.028 

0.025 

LER 

0 

0 

8  exhibiting  maximum  refinement  at  the  wing  tip 
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5.  RESULTS 

The  results  of  this  STTR  Phase  I  are  organized  as  follows.  Convergence  results  are  first  presented, 
followed  by  computational  results.  The  computational  results  include  comparisons  between  methods 
based  on  single  fidelity  level  results  (Section  5.2.1),  application  of  multifidelity  optimization  to  model 
problems  (Section  5.2.2),  and  multifidelity  multidisciplinary  design  problems  involving  aeroelastic  wings. 

5.1  RIGOROUS  CONVERGENCE  RESULTS 

This  section  describes  rigorous  convergence  of  the  methods  that  we  call  “hybrid”  methods  (see 
Section  4.3.5).  The  convergence  proof  is  obvious  and  builds  on  the  multifidelity  convergence  results 
from  Rodriguez,  Renaud,  and  Watson  (1998).  Those  authors  described  rigorous  results  and  showed 
that  their  trust  region  method  converged  under  the  assumptions  they  defined.  Their  results  hold  for 
deterministic  optimization  only  and  are  typical  perhaps  of  nonlinear  programming  convergence  to 
local  minima  results.  In  general,  with  only  a  few  conditions  any  downhill  search  converges.  The 
contribution  here  relates  to  the  multifidelity  aspect  of  the  modeling. 

After  clarifying  their  conditions  and  the  trivial  convergence  of  hybrid  methods  we  discuss  how  some  of 
the  Rodriguez,  Renaud,  and  Watson  (1998)  conditions  might  be  relaxed  and  computationally  efficient 
methods  might  be  developed  with  rigorous  guarantees.  The  hybrid  methods  are  assumed  to  be  likely  to 
offer  global  results  similar  to  other  methods  based  on  global  methods  and  convergence  similar  to  trust 
region  methods.  Yet,  it  is  believed  further  method  development  will  be  needed  to  achieve  both  global  and 
rigorous  convergence  results  and  the  computational  efficiency  of  nonhybrid  methods. 

5.1.1  CONVERGENCE  CONDITIONS 

Trust  region  augmented  Lagrangian  multifidelity  methods  have  been  proven  to  converge  to  local 
optima  of  the  highest  fidelity  system  (Rodriguez,  Renaud,  and  Watson,  1998).  Yet,  their  assumptions 
include  that  the  systems  are  deterministic  and  that  the  fidelity  level,  t| /,  can  be  continuously  adjusted. 

In  the  context  of  multiple  local  minima  they  can  be  considered  inferior  since  sequential  optimization 
methods  using  global  metamodels  have  been  demonstrated  on  test  problems  to  converge  efficiently  to 
the  global  optimum. 

The  purpose  of  this  section  is  to  clarify  the  conditions  used  to  prove  the  trust  region  convergence.  As  a 
result,  reviewing  the  assumptions  from  Rodriguez,  Renaud,  and  Watson  (1998)  clarifies  the  convergence 
criteria  for  hybrid  methods.  It  also  illuminates  the  conditions  and  issues  for  future,  more  efficient  hybrid 
optimization  method  convergence.  Because  the  notation  used  is  somewhat  different  than  in  preceding 
sections  and  because  of  the  length  of  the  material,  the  list  of  underlying  assumptions  is  relegated  to 
Appendix  B. 

5.1.2  CONVERGENCE  RESULTS 

Clearly,  if  a  method  can  be  applied  starting  at  any  point  and  converge,  convergence  is  guaranteed  starting 
from  the  output  of  a  SRBO  or  MFSRBO  method.  Therefore,  simple  hybrid  methods  converge  if  the 
method  used  in  the  last  phase  converges. 

Nontrivial  combinations  of  sequential  radial  basis  optimization  (SRBO)  and  trust  region  components  can 
certainly  be  developed  that  offer  computational  efficiency  advantages  and  proven  convergence  benefits. 
The  key  issues  are  the  issues  for  virtually  all  nonlinear  programming  methods:  (1 )  maintaining  of  positive 
probability  of  improvement  and  (2)  converging  only  when  the  conditions  for  a  local  minimum  (KKT)  are 
met.  Also,  the  proofs  in  Rodriguez,  Renaud,  and  Watson  (1998)  are  presented  in  terms  of  a  series  of 
paired  levels  of  fidelity.  It  seems  likely  that  the  key  steps  can  be  revised  so  that  the  fidelity  levels  can  be 
discrete  and  the  proof  still  applies. 
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5.2  COMPUTATIONAL  RESULTS 

The  MFSRBO  algorithm  outlined  above  was  developed,  tested,  and  refined  based  on  a  numerical  test  bed 
consisting  of  eight  different  optimization  problems.  The  main  characteristics  of  these  eight  problems  are 
outlined  in  Table  7. 


Table  7.  Overview  of  numerical  test  cases. 


Multi 

- 

Design 

Number  of  Constraints 

Test  Problem  Name 

Goal  /  Comment 

Disciplinary 

Fidelity 

Variables 

Aero. 

Structural 

Six  hump  Branin  function 

Unconstrainted 
minimization, 
exhaustive  grid- 
based  search, 
comparison  to  other 
methods 

no 

no 

Xi,  x2 

n/a 

Rosen  brock  function  & 
approximation 

Unconstrainted 

multifidelity 

minimization 

no 

✓ 

Xi,  x2 

n/a 

Symmetrized  Rosenbrock 
function  &  approximation 

Multifidelity 
optimization  with  two 
distinct  optima 

no 

✓ 

Xi,  x2 

n/a 

Aeroelastic  Wing  Design 

1 

a0 

1 

1 

“ 

2 

Minimize  induced 

a0,  da/ds 

1 

1 

“ 

3 

drag,  subject  to  lift 
and  stress 
constraints 

✓ 

✓ 

a0,  da/ds, 
Bo 

1 

1 

“ 

4 

a0,  da/ds, 
Bo,  A 

1 

1 

UAV  Aeroelastic  Wing  Design 

Optimize  stability, 
subject  to  trim, 
payload,  and  stress 
constraints 

✓ 

✓ 

a0,  da/ds 

2 

1 

With  the  exception  of  the  first  test  problem  (Section  5.2.1),  the  surrogate  search  used  in  each  of  the 
optimizations  amounts  to  a  multistart  steepest  ascent9  scatter-and-poll  strategy  (Queipo  et  al.,  2005)  for 
calculating  arg[max[EI  (x ,/)))  . 

5.2.1  SINGLE  FIDELITY  LEVEL  RESULTS 

The  single  fidelity  results  presented  in  this  section  are  almost  the  same  for  the  test  functions  the  OSU 
members  of  our  Team  considered  as  the  results  in  Huang,  Allen,  Notz,  and  Miller  (2006)  and  Huang, 
Allen,  Notz,  and  Zheng  (2006).  This  follows  because  in  virtually  all  of  the  test  problems  good 
solutions  were  obtained  from  the  initial  design  portions  which  were  identical  between  the  Kriging 
model  and  radial  basis  function-based  approaches.  In  all  of  the  problems  in  Huang,  Allen,  Notz,  and 
Miller  (2006),  for  example,  fewer  than  25%  of  the  evaluations  were  after  the  initial  design  evaluations 
First,  we  consider  the  single  fidelity  noisy  case  for  simplicity. 

In  the  results  shown  below  the  Sequential  Radial  Basis  Optimization  (SRBO)  method  is  empirically 
compared  with  three  alternative  approaches  from  the  literature.  The  first  alternative  considered  is  the 
Simultaneous  Perturbation  Stochastic  Approximation  (SPSA)  from  Spall  (1998).  The  second 
approach  is  the  Revised  Simplex  Search  (RSS)  procedure,  a  variant  of  the  Nelder-Mead  method, 
proposed  by  Humphrey  and  Wilson  (2000).  As  the  original  RSS  procedure  was  not  designed  to 
handle  constraints,  we  modify  it  so  that  whenever  a  point  is  infeasible,  the  algorithm  selects  the 

9  see  Appendix  A  for  further  detail 
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nearest  feasible  point  instead.  The  third  approach  is  the  DIRECT  method  developed  by  Gablonsky 
(2001)  etal. 

Here,  we  focus  on  a  single  fidelity  system  with  noise.  In  this  case,  the  test  function  is  the  six  hump 
camelback  function  from  Branin  (1972)  with  ~  N(0,  0.122)  noise  added.  The  function  is: 

/ (x)  =  4  xi2  -  2. 1  xi4  +  1/3  xi6  +  xi  X2  -  4  X22  +  4  X24  (18) 

-1.6  <xi<  2.4,  -0.8  <  X2<  1.2 

x*  =(0-089,  -0.713)  and  (-0.089,  0.713),./*  =  -1.03 

In  this  problem,  there  are  six  local  and  two  global  optima.  For  the  single  fidelity,  noisy  case,  the  expected 
improvement  function  (Eq.  (11))  simplifies  with  CCi  =  a3  =  1. 

Figure  3  plots  two  randomly  selected  sequential  radial  basis  optimization  (SRBO)  runs  overlaid  on  the 
results  of  four  alternative  methods.  The  SRBO  method  was  applied  using  y  =  0.5  and  r  =  3,  i.e.,  the  last 
three  runs  are  used  for  the  residual.  The  results  indicate  that  SRBO  methods  converge  quickly  like 
sequential  Kriging  optimization  (SKO)  methods.  There  might  be  an  issue  with  convergence  although  that 
might  relate  to  the  100x100  grid  used  for  expected  improvement  optimization.  It  is  also  possible  that  the 
convergence  is  limited  by  the  constant  y,  which  explains  why  our  Team  varied  y  in  our  aerospace  wing 
case  study  (Section  5.2.3). 


SRBO  (1) 
SRBO  (2) 
SPSA(1) 
SPSA  (2) 
RSS  (1) 
RSS  (2) 
DIRECT  (1) 
DIRECT  (2) 
SKO  (1) 
SKO  (2) 


Figure  3.  Results  from  two  randomly  selected  runs  of  each  method  on  the  six  hump  camel  back  function. 

(Vertical  axis  corresponds  to  the  function  value.  Horizontal  axis  corresponds  to  the  number  of 
function  calls). 
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5.2.2  MULTIFIDELITY  TEST  PROBLEM 

The  application  of  the  augmented  dimensionality  multifidelity  sequential  radial  basis  optimization 
(MFSRBO)  to  Rosenbrock's  function  is  considered  next.  Similarly  to  Eldred  et  al.  (2004),  we 
consider  the  minimization  of  Rosenbrock's  function  (Figure  4): 

f  hf  =  100(x2  —  x?)2  +  (1  —  Xj)2  (19) 


f  HP  is  the  true  (high-fidelity)  function  to  be  minimized.  It's  low-fidelity  surrogate  f  LF  is  defined  as 


fLP  =  100(x2  -  x\  +  0.2)2  +  (0.8  —  Xj)2 


(20) 


Figure  4.  Isocontours  of  the  high-fidelity 
Rosenbrock  function. 
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Figure  5.  Comparison  between  the  Rosenbrock 

function  and  its  low-fidelity  counterpart  in 
the  vicinity  of  their  optima. 


While  the  two  functions  have  similar  structure,  the  high-fidelity  function  has  its  minimum  at  (  x1  —  1  , 
x2=l  X  the  low-fidelity  one  at  (  x1  =  0.8  ,  x2  =  0.44  ).  The  two  functions  are  not,  in  the  strictest 
sense,  weakly  consistent,  but  the  problem  has  favorable  structure.  In  the  immediate  vicinity  of  the 
high-fidelity  optimum  fHF  =  0  but  f ,,,  »  4  and  the  difference  between  the  two  models  exhibits 
strong  gradients,  as  illustrated  in  Figure  5. 


Because  the  methods  described  in  Section  4.2.2  integrate  data  at  all  levels  of  fidelity  to  build  a  global  radial 
basis  function  metamodel,  it  is  possible  to  use  the  associated  Multifidelity  Expected  Improvement  function 
(Eq.  (1 1))  as  the  criterion  determining  the  location  and  fidelity  level  of  subsequent  evaluations: 

Cost 


EI(x,k)  *  E \max [fm{x  )  -  fm{x)  -  am{x)e ,0 U.corr  fk(x),fm(x) 


1  - 


o\ 


Voa(*) 


+  CT 


0.A-  / 


Cost, 


where  fm{x)  is  the  high-fidelity  prediction  at  point  x  ,  and  x*  is  the  current  “effective  best  solution” 
defined  in  Eq.  (14).  &  20  k  is  the  unexplained  variance  of  the  radial  basis  function  model  at  level  k  ,  and 

(f2(x)  is  the  variance  of  f  k{x)  .  The  correlation  term  corrf  k(x)  ,f  ni(  x)  is  used  to  discount  the 
expected  improvement  when  an  evaluation  from  a  lower  fidelity  surrogate  is  used.  The  term  involving  & 2  k 
accounts  for  diminishing  returns  of  additional  replicates  as  the  prediction  becomes  more  accurate,  and  the 
cost  ratio  adjusts  the  sampling  strategy  according  to  the  evaluation  costs  (Huang  et  al.,  2006).  The  location 
and  fidelity  of  the  next  evaluation,  x„+1  and  kn+1  are  given  by  maximizing  EI(x,k)  ,  i.e.: 
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(xn+1,kn+1)  =  arg[max[EI{x  ,l)\)  (21) 

The  effective  best  solution  x*  isdeftnedby  x*  =  arg\max\u{x )))  where  u(x  j  =  —fm(x)  —  /jsm(x) 
expresses  the  willingness  to  trade  one  unit  of  the  predicted  objective  for  p  units  of  the  standard  deviation  of 
the  prediction  uncertainty  (Schenk  et  al.,  2005).  Also,  the  expectation  in  EI(x  ,k)  is  conditional  given  the 
past  data  and  given  estimates  of  the  correlation  parameters.  Thus,  the  expectation  is  computed  by  integrating 
over  the  distribution  of  f  J  x  )  ,  with  fj  x  )  fixed.  Specifically,  assuming  a  normally  distributed  error,  the 
probability  of  improvement  P[fm(x)<  /m(x*)]  is  defined  as 

/«(*') 

P[fm(x)<fm{x*)]  =  — j-r  J  4>(y-fm(x);sm{x))dy  (22) 

—oo 


where  <f)  is  the  probability  density  function: 


(t>(y;a)  = 


1 


V(2 


exp 


TT 


y 


2  cr" 


(23) 


As  pointed  out  in  Forrester  et  al.  (2009),  the  expected  improvement  is  the  first  moment  of  Eq.  (22),  i.e., 


El 


1 


/„(*) 


/„(*)</„(*: 


sjx 


f  (fm(x)-y)-<f>(y-fm(x);sjx))dy 


(24) 


After  integration  by  parts,  the  expectation  can  be  calculated  analytically  as  follows: 

rnax(7m(x  )-/m(x),0)]  =  (f  m{x)-f  Jx)).y{f  m{x)-f  m{x);s  m{x) 

+  s  (x).<f)if m(x*)—f  (x);s  (x) 


(25) 


where  cf>  is  the  nonnal  probability  density  function  defined  above,  and  (//  is  the  normal  cumulative 
distribution  function: 


v(y;<r)  =  \ 


1  +  erf 


y 


v  2<x 


(26) 


For  the  multifidelity  Rosenbrock  function,  the  initial  design  consisted  of  four  points  (the  four  corners  in 
Figure  4),  three  of  which  were  low-fidelity.  The  high-fidelity  corner  was  chosen  randomly.  An  assumed 
data  point  uncertainty  (  cr0  in  Eq.  (8))  of  0.01  was  used.10  A  typical  optimization  result  is  shown  in 
Figure  6  which  depicts  the  search  path  history.  In  Figure  6  the  low-fidelity  function  evaluations  are 
identified  by  the  small  symbols  and  the  high-fidelity  ones  with  the  larger  symbols. 

The  optimization  correctly  identifies  the  high-fidelity  optimum  in  a  relatively  small  number  of  high- 
fidelity  function  evaluations.  For  reference,  Figure  7  shows  the  low-fidelity  subset  of  Figure  6,  super¬ 
imposed  on  isocontours  of  the  low-fidelity  function.  This  picture  helps  understand  the  role  of  the  low- 
fidelity  evaluation  as  an  aid  to  sample  the  design  landscape.  Note  in  particular  the  accumulation  of  points 
near  the  high-fidelity  optimum  (  x1  —  1  ,  x2  =  1  ),  not  the  low-fidelity  one  (  x1  =  0.8  ,  x2  =  0.44  ). 

In  the  results  depicted  in  Figures  6  and  7  the  assumed  cost  ratio  between  the  two  fidelity  levels  was  3/2. 
As  expected,  the  multifidelity  expected  improvement  function  effectively  drives  the  objective  towards  the 
high-fidelity  optimum  and  takes  advantage  of  the  presumed-low-cost,  low-fidelity  evaluations  to  find  a 
positive  direction  of  improvement.  The  approach  is  successful  because  the  cumulative  radial  basis 
function  metamodel  is  able  to  “learn”  the  relationship  between  the  low-  and  high-fidelity  projections.  In 
an  ideal  sense,  this  nonlinear  relationship  must  be  accurately  defined  using  as  few  (expensive)  high- 

10  for  purposes  of  comparison  with  the  Branin  function  example  considered  in  the  previous  section,  this  data  point 
uncertainty  can  be  thought  of  as  ~  N(0,  0.01)  “virtual”  noise 
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fidelity  observations  as  possible.  Our  empirical  observations  suggest  that,  because  the  uncertainty  on  this 
relationship  remains  large  if  the  number  of  high-fidelity  observations  is  too  few,  there  is  a  practical 
tradeoff  which  must  be  considered  in  the  initial  design.  For  the  majority  of  our  experiments,  the  ratio  of 
high-fidelity  to  low-fidelity  initial  count  was  around  25%. 
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Figure  6.  Multifidelity  search  path  and  isocontours 
of  the  high-fidelity  Rosenbrock  function. 


Figure  7.  Low-fidelity  subset  of  multifidelity  search 
path  with  isocontours  of  the  low-fidelity 
Rosenbrock  function. 


To  further  demonstrate  the  method  on  a  problem  exhibiting  more  than  one  local  optimum,  we  considered 
the  following  dual  optimum  problem,  constructed  by  smooth-symmetrizing  the  Rosenbrock  function 
about  the  origin  as  follows: 


h(x1,x2)  =  ri.f{xux2)  +  (1-/7 )./(— Xj,— x2) 
p  =  ^(tanh[20(x1+x2)]+l) 


(27) 


In  Eq.  (27)  /  designates  either  / HF  or  f  LF  as  appropriate. 


As  in  the  previous  example,  the  initial  design  consists  of  the  four  corner  points  (  x1  —  ±  2  ,  x2  —  ±2  ), 
only  one  of  which  is  evaluated  at  the  high-fidelity  level. 

Figure  8  shows  the  multifidelity  search  path  obtained  by  MFSRBO  with  an  assumed  3/2  cost  ratio  and  an 
assumed  cr0  =  0.01  .  As  in  Figure  6,  the  small  symbols  in  the  figure  represent  the  results  of  low- fidelity 
analyses,  and  the  large  symbols  correspond  to  high-fidelity  analyses.  It  is  noteworthy  that,  regardless  of  the 
fidelity  at  which  they  are  calculated,  the  accumulated  designs  find  both  optima  of  the  high-fidelity  function, 

h  HF  ■ 


The  convergence  history  is  illustrated  in  Figure  9,  which  depicts  in  semi-log  scale  the  values  hHF{x  1 ,  x2) 
of  the  high-fidelity  infills  in  order  of  appearance.  In  order  for  the  optimization  to  be  able  to  both  explore  the 
design  landscape  and  exploit  the  narrow  minima,  it  is  advantageous  to  relax  the  metamodel  stiffness  (i.e., 
increase  y  (Section  4.2.2. 1)  or  y,  (Section  4.2.2.2))  as  the  design  progresses.  It  its  present 
implementation,  this  adaptivity  is  based  on  a  simple,  heuristic  algorithm.  This  algorithm  increases  the  y, 
parameters  by  a  small  factor 


NEW 

Yi 


(1+T, 


OLD 


Yi 


(28) 
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anytime  the  rank  of  the  linear  model  design  matrix  X  recedes  below  a  given  percentage  of  its  full  rank: 

22^  s  T,  ,29) 

For  example,  the  results  of  Figures  8  and  9  were  obtained  using  Tj  =  0.11  and  t2  =  0.94  . 


Figure  8.  Multifidelity  search  path  superimposed 
with  isocontours  of  the  high-fidelity 
symmetrized  Rosenbrock  function. 


Figure  9.  History  of  high-fidelity  infills  corres¬ 
ponding  to  the  mulitifidelity  optimization 
search  path  shown  in  Figure  8. 


5.2.3  AEROELASTIC  WING  DESIGN  PROBLEM 

Two  aeroelastic  wing  design  problems  were  considered  as  part  of  this  investigation.  The  first,  designated 
“baseline  problem,”  consisted  of  an  aspect  ratio  10  unswept  rectangular  or  slightly  tapered  wing  designed 
for  minimal  drag  under  load.  The  second  problem,  referred  to  as  the  “UAV  wing  design,”  was  aimed  at 
maximizing  stability  and  was  characterized  by  the  incorporation  of  multiple  constraints  and  a  noisy 
objective  function. 

5.2.3. 1  Baseline  problem 

A  baseline  problem  was  defined  for  the  purpose  of  developing  and  demonstrating  the  MFSRBO  approach 
and  its  ability  to  represent  and  integrate  information  from  each  discipline  and  model.  Computational  fluid 
dynamics  (CFD)  and  computational  structural  mechanics  (CSM)  were  initially  targeted  as  examples  of 
high-fidelity  models  for  the  two  interacting  disciplines  (fluids  and  structures).  However,  the  practical 
realities  of  the  Phase  I  time  frame  and  budget  focused  our  attention  instead  on  the  following  analysis 
methods  and  fidelity  levels  (see  Section  4.4): 


Table  8.  Overview  of  computational  methods  used. 


MODELS 

Fluids 

Structures 

High 

Fidelity 

OPTMIS/MISDL  (NEAR),  intermediate-level  aerodynamic  prediction  code 
based  on  panel  methods  and  other  classical  singularity  methods  enhanced 
with  models  for  nonlinear  vortical  effects.  Applicable  to  subsonic  and 
supersonic  flight  vehicles  including  aircraft,  UAVs,  missiles,  and  rockets. 

CNEVAL  FEM  method; 
provides  weights, 
displacements,  stresses,  and 
modal  frequencies 

Low 

Fidelity 

Lifting  line  theory,  supplemented  with  induced  drag  calculation. 

Low  degree-of-freedom 
version  of  the  CNEVAL  FEM 
model. 
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An  analysis  of  baseline  problem  candidates  led  to  the  selection  of  the  minimum  drag  design  of  a  wing  subject 
to  minimum  lift  and  static  aeroelastic  constraints  as  the  focus  of  the  Phase  I  demonstration.  This  benchmark 
problem  is  derived  from  Robinson  et  al.  (2006).  It  seeks  to  minimize  the  induced  drag  C D  by  changing  the 
pre-load  twist  distribution  (parameterized  here  by  the  root  chord  angle  of  attack  «0  and  its  spanwise  first  and 
second  derivatives,  daldy  and  d2aldy2)-  Two  additional  variables  are  introduced:  the  root  chord 
maximum  structural  thickness  90  =  TH  (y  =  0)  ,  and  the  taper  ratio  A  =  cljplcroot  of  the  unswept  flexible 
wing.  The  wing  span  is  maintained  constant,  y max  =  5  ft ,  and  the  wing  area  is  kept  fixed  at  S  =  5  ft2  ■ 

The  induced  drag  was  minimized,  subject  to  minimum  lift  and  maximum  stress  constraints  as  follows: 


minimize 

{a0,daldy,90,A} 

subject  to 


C 


D,i 


C  r  ^  payload  +  Wwjnp)l(  S  rpf) 


wing ’ 


ref  ’ 


°Vm  /  CrVM,max<  ^ 


(30) 


where  crVM  is  the  Von  Mises  stress,  W  is  weight,  qx  is  dynamic  pressure,  and  Sref  is  the  wing  reference 
area.  For  this  investigation,  the  freestream  angle  of  attack  is  zero,  the  freestream  Mach  number  is 

=  0.2  ,  the  dynamic  pressure  is  q.,  =  40.76  psf  ,  and  the  assumed  payload  is  W  payload  =  75  lbs  . 


Figure  10.  Low-  and  high  -fidelity  static  aeroelastic  deformation  prediction.  (Top:  aerodynamic 
loading  (planform  view).  Bottom:  streamwise  view  of  combined  wing  twist  and 
bending  (not  to  scale)). 


A  typical  comparison  of  low-  and  high-fidelity  analyses11  is  given  in  Figure  10, 12  corresponding  to  oc0  =  7° , 
daldy  =  -0.8"/A'  ,  d2  al dy2  =  0°  ff 2  ,  0n  =  0.045  ft ,  and  A  =  0.8  using  the  material  and  cross- 
sectional  properties  shown  in  Section  4.4. 

One-,  two-,  three-,  and  four-dimensional  variations  on  the  minimum  drag  aeroelastic  wing  design  problem 
were  investigated  as  part  of  this  effort.  To  better  understand  how  the  optimization  method  performed  on 
these  cases,  let  us  consider  the  two-dimensional  results  (Table  7,  subcase  2)  in  which  a  fixed  root  chord 
maximum  thickness  90  =  0.04  ft  and  fixed  taper  ratio  A  =  1  are  used. 


11  The  cost  ratio  between  high-  and  low- fidelity  analyses  was  5:1 

12  “a”  in  Figure  10  denotes  the  ratio  aVMlcrVM  max 
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Figure  11.  Sample  multifidelity  search  path  for  two- 
dimensional  aeroelastic  wing  design 
problem. 


Figure  12.  A  different  instantiation  of  the  problem 
solved  in  Figure  11  (different  initial 
conditions). 


Figures  11  and  12  provide  two  examples  of  the  MFSRBO  search  path  when  the  optimization  is  initiated 
from  two  randomly  chosen  Latin  Hypercube  design  of  experiments,  each  including  one  high-fidelity  and 
three  low-fidelity  solutions  (see  Table  9). 


Table  9.  Listing  of  initial  conditions 


Point 

Initial  Design  Used  in  Figure  11 

Initial  Design  Used  in  Figure  12 

Number 

«o(°) 

dalds  (°) 

Fidelity  level 

«o(°) 

dalds  (°) 

Fidelity  level 

1 

8.405 

-2.453 

Low 

3.021 

-9.752 

High 

2 

0.784 

-4.876 

Low 

7.487 

-1.617 

Low 

3 

4.690 

-6.679 

High 

9.019 

-6.821 

Low 

4 

5.851 

-8.127 

Low 

0.934 

-3.438 

Low 

The  following  computational  analysis  uncertainties  (  <J0  in  Eq.  (8))  were  used  uniformly  for  the  OPTMIS 
outputs: 


Table  10.  OPTIMS  output  uncertainties 


Cl 

Cdj 

Wwing  (lb) 

CTVM/CTVM,max 

(To 

0.0002 

0.002 

0.03 

0.02 

The  results  of  Figures  11  and  12  are  shown  as  a  function  of  the  two  design  variables  a0  and  d  aids  ,  where 
s  is  the  nondimensional  span  variable,  s  =  yl  ymax  .  The  background  contour  lines  represent  the  high- 
fidelity  primary  (unconstrained)  objective,  C D  , .  The  large  symbols  denote  high-fidelity  analyses,  while  the 
small  symbols  represent  the  low-fidelity  analyses.  Each  symbol  is  colored  based  on  the  objective  function, 
which  is  the  induced  drag  C  D  l ,  modified  with  penalty  functions  for  violating  lift  and/or  stress  constraints 
(themselves  indicated  by  dashed  lines).  Specifically,  the  constraints  are  incorporated  into  the  objective 
function  using  the  following  penalty  formulation  (using  the  previous  values  for  b/payload  ,  and  Srej  ): 
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/  =  C^.+  10 


max 


75  + 

-  Cr  ,0 

+  20 .  max 

40.76x5 

L  3 

cr 


VM 


cr 


-  0.5,0 


VM ,  max 


(31) 


where  W  wing  is  the  structural  weight  of  the  wing,  expressed  in  pounds.  The  process  of  maximizing  the 
multifidelity  expected  improvement  function  (11)  involves  calculating  (15),  which  requires  keeping  track  of 
the  current  best  effective  solution  x  =  [a0  ,(d cx/dy  )  ,00,A  ]  ,  and  the  ability  to  estimate  the  objective 
function  /  and  its  variance  s2  at  any  design  point  x  for  any  level  of  fidelity. 


A  more  elegant  way  of  handling  the  constraints  in  a  consistent  and  fully  probabilistic  manner  is  described  in 
Forrester  et  al.  (2009).  The  detailed  implementation  of  this  method  applied  to  the  present  problem  is 
described  in  Appendix  C. 


Figure  13.  Quantification  of  multifidelity  vs.  high-fidelity-only  optimization  performance  via  aggregated 
cumulative  distribution  functions.  The  primary  objective,  Cdj  is  shown  on  the  horizontal 
axis.  The  vertical  axis  is  the  numerically  computed  probability  of  finding  feasible  high-fidelity 
solutions  exceeding  the  target  on  the  x  axis. 

Statistical  Results 

It  is  important  to  realize  that,  due  to  the  inherent  coupling  between  the  optimization  path  and  the  sequential 
radial  basis  function  metamodel  being  created  at  the  various  stages  of  the  optimization,  the  results  depend  on 
the  initial  conditions  of  the  optimization  (this  is  true  whether  using  multifidelity  or  single  fidelity 
optimization).  Figure  12,  for  example,  presents  a  different  data  profile  for  the  same  nominal  problem  as 
Figure  11,  the  only  difference  being  the  initial  conditions. 
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For  this  reason,  when  comparing  perfonnance  metrics,  it  is  particularly  important  to  expand  results  such  as 
Figures  11  and  12  and  to  consider  ensembles  in  a  statistical  sense  in  order  to  draw  meaningful  conclusions. 
Thus,  using  multiple  realizations  with  random  initial  conditions  based  on  a  Latin  Flypercube  design  of 
experiments,  optimization  results  were  aggregated  and  analyzed  as  a  group.  When  comparing  multifidelity 
vs.  high-fidelity-only  optimization,  care  must  also  be  exercised  so  that  the  exact  same  design  variables  are 
used  in  the  initial  design. 

The  details  of  this  aggregate  analysis  are  as  follows.  The  candidate  designs  harvested  as  a  result  of  each 
optimization  are  first  filtered  so  that  only  the  high-fidelity  (HF)  solutions  that  do  not  violate  the  constraints 
are  considered. 13  These  HF  nonviolators  are  then  sorted  based  on  the  primary  objective  (low  CDj  in  the 
present  case).  Aggregate  perfonnance  profiles  collected  in  this  tnamrer  can  then  be  analyzed  in  tenns  of  the 
induced  drag  cumulative  distribution  function  (CDF)  of  the  wing  designs.  The  example  shown  in  Figure  13 
corresponds  to  8  realizations  of  the  three-dimensional  optimization  (Table  7  subcase  3).  The  stopping 
criterion  was  based  either  on  failure  to  improve  the  El  function  (Eq.  (11)),  or  reaching  an  imposed  practical 
maximum  of  24  high-fidelity  solutions.  The  optimal  designs  for  this  case  are  given  in  the  table  below: 


Method 

«o(°) 

daldy  (° /L1) 

00  ( ft ) 

CL 

cDJ 

Wwing  (lbs) 

^VM 

MFSRBO 

6.624 

-0.704 

0.0392 

0.4624 

6.484  x  lO'3 

18.87 

0.442 

HF-only 

6.418 

-0.603 

0.0384 

0.4643 

6.542  x  lO'3 

18.64 

0.471 

Figure  13  shows  that,  by  using  multifidelity  optimization,  there  is,  for  example,  a  37%  probability  of 
achieving  a  drag  coefficient  less  than  0.01,  versus  less  than  4%  with  high-fidelity  optimization  only.  Thus, 
for  a  given  computational  budget,  there  is  an  order  of  magnitude  greater  probability  of  finding/approaching 

the  high-fidelity  optimum  bv  using  multifidelitv  optimization. 14 


o 
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Figure  14.  Relative  efficiency  between 

multifidelity  and  high-fidelity-only 
optimization  as  a  function  of  design 
space  dimensionality. 


Neglecting,  for  the  sake  of  simplicity,  the  cost  of  the  low- 
fidelity  analyses  and  the  cost  of  the  optimization 
calculations  per  se,  the  average  ratio  of  the  cumulative 
distribution  functions  provides  an  indication  of  the  relative 
performance  of  these  methods.  Let  us  define  the  CDF  area 
ratio  as  the  ratio  between  the  areas  under  the  CDF  curves 
(Figure  13)  for  CDJ  <  0.01  . 

Figure  14  shows  the  CDF  area  ratio  as  a  function  of  the 
search  space  dimensionality  of  the  problem.  The  three-  and 
four-dimensional  cases  are  believed  to  be  more  typical  of 
the  expected  performance  ratio.  Our  observations  suggest 
that  the  rcplicatc-dctcrrcnt  term  (13),  may  have  had  an 
important  effect  in  the  lower  dimensional  experiments,  due 
to  the  effectively  higher  point  density  in  those  cases.  The 
one-dimensional  case  was  obtained  by  varying  a0  only. 
(The  so-called  1. 5-dimensional  case  represents  a  case 
where  the  constraints  were  such  that  only  a  narrow  sliver  of 
feasible  solutions  existed  in  the  2-D  space,  thus  making  the 
space  effectively  quasi-one-dimensional).  While  limited  in 


scope  and  strictly  empirical,  the  results  of  Figure  14  suggest  the  possibility  that  the  MFSRBO  method  may 
be  increasingly  beneficial  as  the  dimensionality  of  the  problem  increases,  although  this  question  will  be  left 
for  future  research. 


13  In  this  study,  this  was  typically  repeated  for  8  to  16  realizations. 

14  This  is  consistent  with  the  prior  MFSKO  experience  of  the  OSU  members  of  our  Team,  which  suggests 
dramatic  benefits  when  the  cost  differential  is  high  and  the  systematic  errors  are  small. 
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5.23.2  UAV  wing  design 


The  baseline  problem  described  above  was  extended  to  the  case  of  an  aeroelastic  swept  (UAV-type)  wing  at 
=  0.2  .  The  wing  planform  is  a  modified  version  of  the  European  SACCON  stability  and  control 
vehicle15  with  identical  span  and  wing  area.  As  in  Table  7  subcase  2  above,  the  design  variables  are  the  twist 
distribution,  characterized  by  tx0  and  da  Ids  .  However,  this  time,  the  design  objective  is  to  maximize  the 
pitch  stability  —dCjd  a  ,  subject  to  minimum  lift,  maximum  stress,  and  angular  trim  constraints  as 
follows: 


minimize 
{ aQ,d  aldy } 

subject  to 


dCJ  da 


CL  ~  payload+WWiney( do 


aVM  t  aVM ,  max 

I  C  I  <0.002 


wing ’ 

<  0.55 


5 


ref' 


(32) 


where  Cm  is  the  pitching  moment  about  the  root  half  chord  location. 


An  example  comparison  of  low-  and  high-fidelity  results  obtained  using  OPTMIS  is  shown  in  Figure  15, 16 
corresponding  to  cxit  =  1 8.5"  ,  daldy  =-4.3° /T1  and  W payload  =75  lbs  for  an  isotropic  material  with 
self-similar  cross-section  (see  Section  4.4)  and  a  spanwise-uniform  maximum  thickness  of  dQ  =  0.025  ft . 


Figure  15.  Low-  and  high-fidelity  static  aeroelastic  deformation  of  UAV  wing. 

(Left:  planform  view  of  loading.  Right:  perspective  view  of  wing  twist 
and  bending  (not  to  scale)). 

The  optimization  methods  and  postprocessing  metrics  were  identical  to  those  used  in  Section  5.2.3. 1 . 
Figure  16  provides  an  illustration  of  the  multifidelity  sequential  radial  basis  optimization  search  path  for  a 


15  NATO  RTO,  AVT- 161  Assessment  of  Stability  and  Control  Prediction  Methods  for  NATO  Air  &  Sea  Vehicles. 

16  “a”  in  Figure  15  denotes  the  ratio  aVMlaVM  max 
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given  instance  of  a  Latin  Hypercube  design  of  experiments  including  one  high-fidelity  and  three  low-fidelity 
solutions  as  the  initial  conditions.  The  flooded  contours  represent  the  primary  objective,  dC  J  d  a  ,  which  is 
noisy,  due  in  part  to  finite  differencing  and  to  Cm  being  itself  a  sensitive  quantity.  As  before,  the  large 
symbols  denote  high-fidelity  analyses,  while  small  symbols  are  used  for  low  fidelity.  Each  symbol  is  colored 
based  on  the  objective  function,  which  is  the  stability  derivative  dC ml da  ,  augmented  with  penalty  functions 
for  violating  either  the  lift,  trim,  or  stress  constraints.  Specifically,  the  constraints  are  incorporated  in  the 
objective  function  using  the  following  penalty  fonnulation  (using  the  previous  values  for  f^payi0ad  >  d,.  > 
and  Sr£f  =  4.15/t2  ): 


/ 


dC, 

d  a 


+  10. max jCJ  —  0.002, 0 


+  10.  max 

+  20.  max 


75  + IT. 


wing 


40.76x4.15 


~CL’0 


VM 


a 


-0.5,0 


VM ,  max 


(33) 


where  TH  win„  is  the  structural  weight  of  the  wing,  expressed  in  pounds. 
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Figure  16.  Sample  multifidelity  search  path  for  UAV  aeroelastic  wing  design  problem. 

Finally,  Figure  17  compares  the  aggregate  cumulative  distribution  functions  for  the  multifidelity  and  single 
fidelity  optimizations,  each  based  on  eight  random  realizations.  The  resulting  CDF  area  ratio,  defined  as  the 
ratio  between  the  areas  under  the  CDF  curves  (Figure  17)  for  dCj da<  —0.0027  ,  suggests  4: 1  odds  in 
favor  of  multifidelity  optimization,  versus  high-fidelity  only  optimization. 
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Figure  17.  Quantification  of  multifidelity  vs.  high-fidelity-only  optimization  performance  via  aggregated 
cumulative  distribution  functions.  The  primary  objective,  dCm/da  is  shown  on  the  horizontal 
axis.  The  vertical  axis  is  the  numerically  computed  probability  of  finding  feasible  high-fidelity 
solutions  exceeding  the  target  on  the  x  axis. 


6.  DISCUSSION 

This  report  has  explored  the  use  of  radial  basis  functions  instead  of  Kriging  models  in  sequential 
multifidelity  optimization.  In  adapting  the  methods,  we  have  addressed  at  least  three  types  of 
unavoidable  challenges  relating  prediction,  error  estimation,  and  the  expected  improvement  function. 
This  investigation  explored  proposing  two  types  of  metamodeling  approaches  which  are  cascading  (the 
most  direct  extension  of  Kriging-based  methods  with  convergent  prediction)  and  augmented 
dimensionality  (which  offers  advantages  for  data  starved  cases). 

Perhaps  the  greatest  challenge  was  developing  a  comprehensive  estimate  of  uncertainty  corresponding  to 
the  approximate  Bayesian  intervals  for  Kriging  functions.  They  are  approximate  because  they  are  based 
on  parameter  estimates.  Here,  four  types  of  approaches  were  considered  with  numerical  results  based  on 
assumed  errors  and  simple  cross-validation.  Finally,  two  alternative  formulations  for  the  expected 
improvement  function  were  described  because  of  the  computational  challenge  of  estimating  correlation 
for  the  radial  basis  function  case. 
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7.  CONCLUSIONS  AND  FUTURE  WORK 

The  preliminary  numerical  results  indicate  that  methods  based  on  radial  basis  functions  can  compete  with 
Kriging  methods  in  relation  to  efficiently  and  repeatedly  identifying  global  optimal  solutions.  As  a  result, 
the  radial  basis  functions  might  seem  to  be  preferred  for  the  following  reasons: 

•  Reduced  computational  overhead  -  Radial  basis  function  estimation  does  not  require  likelihood 
minimization  since  we  have  a  closed  form  for  the  coefficient  estimates.  Thorough  and  accurate 
likelihood  optimization  is  a  difficult  computational  challenge.  The  efficiency  increase  through 
eliminating  one  of  the  two  optimizations  per  iteration  (the  other  is  expected  improvement 
optimization)  can  enable  many  new  applications  of  related  methods.  These  included  so-called 
“on-line”  or  automatic  optimal  decision-making.  For  example,  data  can  stream  in  from  sensors  or 
simulations  running  in  parallel  and  be  quickly  integrated  for  determining  optimal  settings  and 
requests  for  additional  inputs. 

•  Improved  numerical  stability  -  In  the  prior  computational  work  of  the  OSU  members  of  our 
Team,  generated  by  Huang,  Allen,  Notz,  and  Miller  (2006),  we  encountered  several  crashes  from 
numerical  inversions  in  deriving  numerical  results.  These  problems  likely  were  not  fully 
understood  by  us.  However,  our  preliminary  conclusion  was  that  the  errors  occurred  because  of 
the  spacing  off  design  points  near  each  other  and  on  a  lattice  and  could  not  be  avoided  through 
the  use  of  singular  value  decomposition-based  inversion.  By  contrast,  the  radial  basis  function 
fitting  is  based  on  linear  model  estimation  which  is  very  stable  using  singular  value 
decomposition  (Press,  Teukolsky,  Vetterling,  and  Flannery,  2007). 

•  No  stationarity  assumption  -  In  the  Huang,  Allen,  Notz,  and  Miller  method  there  was  a  check 
for  violation  of  the  assumption  that  systematic  errors  have  zero  expected  mean.  The  methods 
included  a  potentially  expensive  checking  and  adjustment  process  which  was  not  fully  reflected 
in  the  measured  computational  performance.  The  methods  based  on  radial  basis  functions  seem 
to  require  no  comparable  steps  which  could  offer  computational  advantages. 

Table  11  qualitatively  compares  alternative  methods  in  relation  to  modeling  approach,  computational 
overhead,  available  convergence  results,  and  test  performance.  The  tentative  conclusion  is  that  further 
hybrid  methods  based  on  radial  basis  functions  are  a  promising  topic  for  research.  The  developed 
methods  will  likely  offer  reduced  modeling  overhead,  theoretical  convergence  results,  and  good  test 
performance  on  numerical  examples  and  real-world  problems. 


Table  11.  Overview  of  alternative  methods.  *Developed  in  this  project..  **Proposed  for  future  development. 


Method 

Modeling 

Modeling  Overhead 

Convergence 

Test  Performance 

SKO  (Jones  et  al.  &  Huang,  et  al.  #1) 

Kriging 

Likelihood  Optimization 

Not-established 

Generally  Global 

SRBO* 

RBF 

LSE 

Not-established 

Competitive 

MFSKO  (Huang.  Allen  et  al.  #2) 

Kriging 

Likelihood  Optimization 

Not-established 

Generally  Global 

Trust  Region  (Rodriguez,  Renaud,...) 

Polynomial 

LSE 

Established 

Not  Competitive 

MFSRBO:  Assumed  Variance* 

RBF 

LSE 

Not-established 

Competitive 

MFSRBO:  Cascading  model* 

RBF 

LSE 

Not-established 

Not  studied 

MFSRBO:  PRESS* 

RBF 

PRESS  Cross  Validation 

Not-established 

Not  studied 

Simple  Hybrid* 

RBF 

LSE 

Established 

Not  Competitive 

Hypothetical  Future  Hybrid  Method** 

RBF 

LSE 

Promising 

Competitive 

The  research  has  also  uncovered  many  remaining  challenges.  These  include: 

•  Developing  globally  convergent,  efficient,  hybrid  methods  -  The  simple  hybrid  methods 
considered  always  converge  to  local  minima  and  might  likely  converge  to  global  minima  under 
the  assumptions  considered.  Yet,  their  computational  performance  has  not  been  studied  and  is 
likely  not  competitive.  Also,  the  assumptions  of  deterministic  responses  and  continuously 
adjustable  fidelity  levels  might  need  to  be  relaxed.  It  is  possible  that  new  classes  of  hybrid 
methods  will  offer  competitive  computational  performance  and  proven  convergence. 
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•  Error  estimation  including  bias  -  The  errors  used  in  the  proposed  MFSRBO  explored  in  our 
computational  studies  are  based  on  regression  variance  errors  only.  This  makes  their 
interpretation  difficult  in  deterministic  cases,  i.e.,  cases  with  perfect  repeatability  error  or  o  =  0. 

In  such  cases,  the  intervals  do  not  directly  correspond  to  any  interpretable  physical  quantity  in  the 
corresponding  problem  and  performance  of  the  methods  likely  is  degraded.  Developing  new 
methods  based  on  measures  of  bias  such  as  regression  Cp  or  diagnostics  proposed  in  Tseng  (2007) 
is  a  promising  avenue  of  research. 

•  Measuring  computational  overhead  -  It  is  clear  that  linear  model  estimation  is  computationally 
superior  to  likelihood  estimation  in  Kriging  models.  Yet,  so  far  we  have  not  measured  the 
computational  advantages.  Also,  we  have  not  tuned  the  methods  to  the  point  at  which  the 
overhead  is  minimal  and  comparable  to  standard  nonlinear  programming  methods.  The 
promising  results  from  the  radial  basis  function  methods  here  suggest  that  such  performance  is 
achievable. 

•  More  thorough  computational  comparisons  -  There  are  several  cases  in  which  further 
comparisons  can  shed  light  on  the  conditions  under  which  radial  basis  function  methods  are 
competitive  or  superior  to  methods  based  on  Kriging  models  and/or  trust  region  methods.  Further 
comparisons  are  possible,  particularly  including  discrete  or  categorical  factors. 

•  Parallelizing  the  methods  -  One  promising  area  that  we  considered  was  the  possibility  to 
parallelize  function  evaluations.  For  example,  it  is  advantageous  to  launch  one  slow,  high- 
fidelity  simulation  simultaneous  with  several  low-fidelity  simulations.  The  structure  of  related 
methods  lends  itself  nicely  to  parallel  processing  with  its  significant  computational  advantages. 
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AFRL 

AR 

arg 

c 

C 

CD,t 

CDF 

CDL 

CFD 

CL 

Cm 

CSD 

CSM 

d 

DIRECT 

DOE 


Air  Force  Research  Laboratory 
Aspect  Ratio  (total  span  over  chord) 

Argument 
Chord  length 
Cost 

Induced  drag  coeffient 
Cumulative  Distribution  Function 
Cross-Disciplinary  Link 
Computational  Fluid  Dynamics 
Lift  coefficient 
Pitching  moment  coefficient 
Computational  Structural  Dynamics 
Computational  Structural  Mechanics 
Design  space  dimensionality 

Optimization  method  developed  by  Gablonsky  et  al.  (2001) 
Design  of  Experiments 
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EGO 

El 

/_ 

Pi 

FEM 

g 

GCV 

HF 

KKT 

KRG 

LER 

LF 

LSE 

MAV 

max 

MDA 

MDA 

MFSKO 

MFSRBO 

min 

MLS 

N[p,c] 

NEAR 

obj 

OSU 

P[] 

Param. 

PRESS 

PRG 

qPararn. 

q°° 

RBF 

ref 

RS 

RSS 

S 

SCO 

SKO 

SPSA 

SRBO 

STTR 

SVD 

SVR 

TH 

“i 

UAV 

Var 

W 

x 

X 

X 


Efficient  Global  Optimization 

Expected  Improvement 

Function 

Force  vector 

Finite  Element  Method 

Constraint 

Generalized  Cross  Validation 
Fligh-Fidelity 
Karush-Kuhn-T  ucker 
Kriging-based  regression 
Leading  Edge  Radius 
Low-Fidelity 
Least  Squares  Estimation 
Micro  Air  Vehicle 
Maximum 

Multidisciplinary  Design  Analysis 
Multidisciplinary  Design  Optimization 
Multifidelity  Sequential  Kriging  Optimization 
Multifidelity  Sequential  Radial  Basis  Optimization 
Minimum 

Moving  Least  Squares 

Random  process  with  normal  distribution  of  mean  p  and  standard  error  a 

Nielsen  Engineering  &  Research 

Objective  function 

Ohio  State  University 

Probability 

Parametric 

Prediction  Error  Sum  of  Squares 

Polynomial  regression 

Quasi  parametric 

Freestream  dynamic  pressure 

Radial  Basis  Function 

Reference  quantify 

Response  Surface 

Revised  Simplex  Search 

Area 

Standard  error  of  function / 

Sequential  Kriging 

Simultaneous  Perturbation  Stochastic  Approximation 
Sequential  Radial  Basis  Optimization 
Small  Business  Technology  Transfer 
Singular  Value  Decomposition 
Support  Vector  Regression 

Maximum  structural  thickness  at  a  given  span  location 

Velocity  vector 

Unmanned  Aerial  Vehicle 

Variance 

Weight 

Streamwise  coordinate 
Design  vector 
Design  matrix 
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y  Spanwise  coordinate 

Y  Dependent  variables  data  vector 

z  Normal  coordinate 

A  Estimated  value 

*  Effective  best  solution 

a  Angle  of  attack 

a0  Root  chord  angle  of  attack 

y  Metamodel  parameter 

T  i  Circulation 

ACp  Differential  pressure  coefficient 

A  Taper  ratio 

p  Fluid  density 

Gvm  Von  Mises  stress 

0O  Root  chord  maximum  structural  thickness 
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APPENDIX  A  -  CORRELATION  COEFFICIENT 

COMPUTATION 


The  objective  function  of  the  optimization  is  the  multifidelity  expected  improvement  function  driver 
EI(x,/)  (see  Eq.  (11)).  At  any  given  point  x  within  the  optimization  search,  the  fidelity-maximization  of 
EI(x,/)  results  primarily  from  the  balance  between  the  cost  ratio  a  3  ( / )  and  the  local  correlation  (/.  {(x  ,1 ) 

between  fidelity  levels  l  and  m,  i.e.,  cc1{x ,  l)  =  corr[/;(x ) ,  f  m  (x )]  .  This  term  is  key  to  the  functioning 
of  the  method:  whenever  lower  and  higher  fidelity  estimates  are  correlated,  the  optimization  favors  the 
lower  cost  function,  and  where  they  are  sufficiently  uncorrelated,17  the  high-fidelity  estimate  takes 
priority. 

In  the  augmented  dimensionality  version  of  our  method,  the  correlation  coefficient  a^x  ,1)  is  computed 

numerically  as  the  correlation  coefficient  between  two  data  streams.  The  two  data  streams  consist  of 
local  sets  {/;( x q) >■■■>  f  iix q)\  and  [f  m(x  o)  >•••  >f  m(x  q)}  ,  where  x0’---’xq\  denotes  a  set  of 
points  in  the  immediate  neighborhood  of  x,  and  /  is  the  metamodel  response.  In  practice,  one  must 
decide  on  what  these  neighbor  points  should  be  (scale  of  coverage,  spatial  pattern,  and  number). 

Our  implementation  “piggy  backs”  on  the  search  stencil  of  the  steepest  ascent,  hill  climbing  algorithm 
used  in  NEAR  RS  (Section  4.3.6).  In  its  simplest  form  (but  not  the  most  efficient),  the  stencil  consists  of 
a  center  point  (  x  0  =  x  )  and  the  vertices  of  a  hypercube  of  size  L,  centered  on  x.  In  our  search  strategy, 
L  is  initiated  at  some  user-prescribed  initial  value  around  the  random  seed  points.18  The  stencil  is  then 
interrogated  at  each  of  the  vertices.  If  an  improvement  is  possible  the  search  location  is  moved  to  the 
vertex  corresponding  to  the  best  improvement  (steepest  ascent),  leaving  L  unchanged.  If,  on  the  other 
hand,  there  is  no  improvement  at  any  of  the  vertices,  then  the  pattern  remains  centered  on  x()  and  L  is 
cut  in  half.  The  algorithm  goes  on  until  reaching  a  maximum  user-prescribed  number  of  recursion  levels 
(effectively  defining  the  maximum  resolution  of  our  search). 

For  the  purposes  of  computing  the  correlation  coefficient,  the  local  length  scale  L  is  a  natural  choice 
for  defining  the  scale  of  coverage  of  the  neighborhood  set  *  o  *  •  •  • >  x  q  J  •  In  this  study,  a 
(/-dimensional  hypercube  of  size  L/2  with  w  =  3  points  per  side  was  used  for  most  of  the  cases, 
providing  reasonably  smooth  estimates  of  the  correlation  coefficient. 


17  (for  a  given  cost  ratio) 

18  The  random  seed  points  are  the  initializations  of  the  multistart  method. 
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APPENDIX  B  -  CONVERGENCE  CONDITIONS 


Trust  region  augmented  Lagrangian  multifidelity  methods  have  been  proven  to  converge  to  local 
optima  of  the  highest  fidelity  system  (Rodriguez,  Renaud,  and  Watson,  1998).  Yet,  their  assumptions 
include  that  the  systems  are  deterministic  and  that  the  fidelity  level,  vp,  can  be  continuously  adjusted. 

In  the  context  of  multiple  local  minima  they  can  be  considered  inferior  since  sequential  optimization 
methods  using  global  metamodels  have  been  demonstrated  on  test  problems  to  converge  efficiently  to 
the  global  optimum. 

The  purpose  of  this  system  is  to  clarify  the  conditions  used  to  prove  the  trust  region  convergence.  As  a 
result,  reviewing  the  assumptions  from  Rodriguez,  Renaud,  and  Watson  (1998)  clarifies  the  convergence 
criteria  for  hybrid  methods.  It  also  illuminates  the  conditions  and  issues  for  future,  more  efficient  hybrid 
optimization  method  convergence.  The  notation  was  somewhat  different  than  in  proceeding  sections.  The 
relationship  is  clarified  below. 


k 

Lagrangian  iteration 

s 

approximate  minimization  iteration 

i,j,  l 

variable  indices 

m 

number  of  inequality  constraints 

n 

number  of  design  variables 

P 

number  of  equality  constraints 

M 

objective  function 

g(x) 

inequality  constraint  vector 

sOO 

y'-th  inequality  constraint 

/?(x) 

equality  constraint  vector 

hi(\) 

/-th  equality  constraint 

rp 

penalty  parameter 

X 

design  vector,  dimension  n 

Xl 

/-th  design  variable 

xu 

upper  bound  vector,  dimension  n 

*U1 

l-th  design  upper  bound 

xL 

lower  bound  vector,  dimension  n 

l-th  design  lower  bound 

B 

approximation  of  the  Hessian 

K 

constraints  residual 

S 

design  space 

P,  Pu  P2,  k 

scalars 

Si,  82 

convergence  tolerance 

To,  Yi,  72,  fi,  p 

trust  region  parameter 

Lagrange  multiplier  vector,  dimension  m+p 
i-th  Lagrange  multiplier 
trust  region  ratio 
fidelity  control 
Euclidean  norm 

Gradient  operator  with  respect  to  design  vector  x 
P[v(x)]  Projection  operator;  projects  the  vector  y  onto  the  set  of  feasible  directions  at  x 

A  trust  region  radius 

Ax  step  size 


X 

h 

P 

V 

V 
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Given  minimization  problems,  authors  made  several  assumptions  as  follow: 


Min:/(x) 
Subject  to 
hi(\)  =  0 ,i  =  1, 
g/(x)  >0,j= 
x£/<x/<xu/,  1  = 


:  equality  const 
:  inequality  const. 
: bounds 


(34) 


(35) 

(36) 

(37) 


(38) 


In  terms  of  these  definitions,  the  nine  key  assumptions  are  as  follows.  It  is  not  clear  whether  all  nine  are 
really  needed  and  they  relate  generally  to  bounds  on  the  roughness  of  relevant  functions.  Notes  follow 
each  assumption.  Overall,  these  assumptions  might  be  considered  relevant  to  many  types  of  modeling 
problems  that  do  not  involve  hits  or  misses  of  moving  components  or  other  forms  of  highly  discontinuous 
physical  behaviors. 

AS1.  Design  space,  S,  has  a  nonempty  interior  and  compact,  i.e.,  it  contains  its  boundary. 

AS2.  Objective  function  (/),  inequality  constraints  (gy),  and  equality  constraints  (hi)  are  twice  continuously 
differentiable  on  an  open  set  containing  the  design  space  S. 

AS3.  The  objective  function  (the  augmented  Lagrangian  function),  ^/i(x),  satisfies  a  Lipschitz  condition. 
In  words,  the  gradients  are  bounded. 

AS4.  The  feasible  set  defined  by  the  inequality  constraints,  equality  constraints,  and  upper  and  lower 
bounds  is  nonempty. 

AS5.  Trust  region  ratio,  ps,  — >  1  as  trust  region  radius,  As  — >  0  or  fidelity  control,  v|/  — >  \| )max. 

Trust  region  “ratio”  is  as  follow  by  the  authors,  5  is  the  iteration  counter  for  the  number  of 
approximate  augmented  Lagrangian  minimizations.  If  ps  <  0,  the  approximation  is  bad  and  the 
trust  region  radius,  As  should  be  reduced.  If /T  >  1,  the  approximation  is  also  bad  but  the  direction 
is  right.  If  psc  (0,  1),  the  radius  can  be  kept  unchanged,  reduced,  or  increased  based  on  how 
close  to  0  or  1  it  is. 


P 


AS6.  Inequality  constraint  vector,  g  and  equality  constraint  vector,  h  satisfy  the  Kuhn-Tucker  constraint 
qualification  or  the  modified  Arrow-Hurwicz-Uzawa  constraints  qualification  at  any  local  minimum. 

•  KKT  constraint  qualification  -  i)  ACQ(Abadie  constraint  qualification)  -  the  tangent  cone  of 
a  feasible  set  at  a  feasible  point  x*  is  the  linearized  cone  at  x*  in  the  tangent  cone,  ii) 
LICQ(Linear  independence  constraint  qualification)  -  the  rows  of  Jacobian  matrices  for 
inequality  and  equality  constraints  are  linearly  independent  (from  a  note  of  The  Karush- 
Kuhn-Tucker  Theorem,  by  Moritz  Kuhn,  2006) 

AS7.  The  Jacobian  matrix  of  the  active  constraints  has  full  rank  at  any  limit  point  x*  of  the  sequence  {xA} 
considered. 
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AS8.ThesetM  —  {{x|0(x)  <  0(x0)]  O  {xL<  x  <  xu} }  with  x°  satisfying  xL,  <  x,  <xUi  (i=  1 is  not 
empty. 

•  It  is  necessary  that  the  solution  of  the  optimization  problem  be  nontrivial. 

AS9.  Approximation  of  Hessian,  ||  B  ||  <  an  assumed  constant. 

•  It  is  necessary  to  assure  that  the  Hessian  (containing  the  second  derivatives)  approximation 
remains  uniformly  bounded. 
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APPENDIX  C  -  DISCUSSION  ON  THE  HANDLING  OF 

CONSTRAINTS 


The  minimum  induced  drag  design  considered  in  the  baseline  problem  of  Section  5.2.3. 1  exercised 
structural  and  aerodynamic  constraints.  Specifically,  the  following  problem  was  solved: 


minimize  C  D  . 

{ cx0  ,daldy,  90,  A } 


subject  to 


CL  ^  ( ' W  payload  +  W wing )  A 


ref ’ 


°Vm  /  <JVM,max<  ^ 


(39) 


The  constraints  were  incorporated  into  the  objective  function  using  a  simple  penalty  method,  specifically 
(  ^payload  =  75  Ibs  ,  =  40.76  psf  and  Sref  -  3  .//  ’  ); 


/  =  C D  j  +  10  .max 


75  +  fK 


wing 


\  40.76x5  j 


-  CT  ,0 


+  20 .  max 


a 


VM 


-  0.5,0 


VM ,  max 


(40) 


where  IT  wing  is  the  structural  weight  of  the  wing,  expressed  in  pounds.  The  process  of  maximizing  the 
multifidelity  expected  improvement  function  (11)  involves  calculating  (15),  which  requires  keeping  track  of 
the  current  best  effective  solution  x*=  \a*0,{d(x/dy)* ,9*0,A*]  ,  and  the  ability  to  estimate  the  objective 
function  /  and  its  variance  s 2  at  any  design  point  x  for  any  level  of  fidelity. 


Each  independent  metamodel  surrogate  for  C D  j  ,  CL  ,  Wwjng  ,  and  (Tvm^(T  vm  .max  yields  an  estimated 
value  (denoted  with  a  ( A )  )  and  standard  error  s  .  The  approach  used  in  this  work  considers  the  primary 
objective  C D  i  as  probabilistic.  However,  for  simplicity,  the  constraints  were  treated  as  deterministic. 
Thus 


/  =  Cn  .  +  10  .max 

J  D  ,i 


15  +  W  ■ 

wing 


40.76  x  5 


+  20 .  max 


a 


VM 


cr 


-  0.5,0 


VM ,  max 


but 


cr 


S(f)  =  s(Cni)  +  0.s(Cr)  +  0.s(W  ■  )  +  0.s( - — ) 


wing’ 


a 


VM ,  max 


(41) 

(42) 


A  better,  fully  probabilistic,  approach  which  takes  into  account  the  probability’  of  constraint  violation  is 
described  in  Forrester  et  al.  (2009),  in  which  the  expected  improvement  is  replaced  by  the  expectation  of 
improving  of  the  objective  f  while  simultaneously  being  a  feasible  design: 


max 


fm(x)-fm{x),0 1  n  {g>  gmin} 


(43) 


In  (43),  g  represents  a  constraint  and  g  >  gmin  represents  the  condition  for  feasible  designs.  In  the  case 
where  primary  objective  and  constraint  are  uncorrelated,  the  expectation  (43)  can  be  expressed  as  a 
simple  product  of  independent  probabilities: 


inaxj7m(x‘)— /m(x),0|  n  {g>  gmin} 


=  E 


max( /J **)-/, „(*).0) 


P  S>  gn 


(44) 
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In  this  fully  probabilistic  approach,  the  penalty-based  objective  function  /  in  Eq.  (40)  is  simply  replaced  by 
/  =  Cdj>  E  max (  /,„(*  )~/m(*)>0)  is  calculated  according  to  (25)  using  sm(x)  =  sm{C D  j;x)  , 


P\gJx)>gmiB\  =  J  <f>(y-gjx);sm(g;x))dy 


=  —  1  +  erf 
2  J 


gjx)-gmin 

,V2.s_(g;x) 


=  V  (gjx)-gwb,;s„(g;x)) 


Thus,  in  the  above  example,  and  provided  the  above-mentioned  independence  assumptions  hold,  (25), 
(41)  and  (42)  are  replaced  by 


£[max(/m(x  )-fpm(x), 0)  n  {g  >gmin}]  = 

[[CD,  ^X^~ED,i  (X)|  ^  —  ( X^’SnSED,i’X )] 

\  m  m  j  \  m  m  / 

+  Stn(C D,i’  x).(t>lc D,i  (X^~ED,i  ( X )»  Sm^D,  i ;X  ))  ] 


\w  75  -lc  (x)  +  ^wing^^ 
j  40.76x5 


>\  Sl(CL’’X)  + 


s2  (W ^m;x) 

m v  wing  ’ 

(40.76  x5)2 


Vvm(x) 


~  0.5 ;s. 


VM.max  /  m 


■;*  ] 


which  explicitly  takes  into  account  the  estimated  variances  sm\CD  ; x )  ,  sm  ,  Sm(W wing’*)  > 

and  sja^la  VM,max  > x )  for  aH  the  dependent  variables  participating  in  either  the  objective  or  the 
constraints. 
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